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A  CONTINUUM  MODEL  FOR  STREAMFLOW  SYNTHESIS 


1.  INTRODUCTION 

Streamflow  evolves  as  a  continuum,  and  is  normally  comprised  of  three 
components:  (1)  surface  runoff,  (2)  interflow,  and  (3)  baseflow.  These 
components  occur  concurrently,  although  their  relative  magnitudes  vary  with  time.  If 
we  consider  a  sudden  burst  of  rainfall,  then  surface  runoff  predominates  during  the 
rising  part  of  the  streamflow  hydrograph,  interflow  predominates  during  the  early  part 
of  its  recession,  and  baseflow  predominates  during  the  delayed  part  of  its  recession. 
The  mechanisms  and,  therefore,  the  governing  equations,  of  these  components  are 
different  but  are  influenced  by  dynamic  interactions  prevailing  between  them. 

Although  streamflow  synthesis  has  long  been  a  subject  of  scientific  inquiry, 
treatment  of  streamflow  as  a  continuum  taking  into  account  dynamic  interactions 
amongst  its  components  has  not  yet  been  fully  developed.  Most  of  the  approaches 
of  streamflow  synthesis  are  based  on  the  concepts  embodied  in  Horton's  infiltration 
theory  of  runoff  (Horton,  1933).  According  to  this  theory,  rainfall  is  absorbed  for 
intensities  not  exceeding  infiltration  capacity,  while  for  excess  rainfall  there  is  a 
constant  rate  of  absorption  as  long  as  the  infiltration  capacity  is  unchanged.  Thus, 
infiltration  divides  rainfall  into  two  parts.  One  part  travels  over  the  surface  giving  rise 
to  surface  runoff,  and  the  other  part  infiltrates  into  the  ground  resulting  in 
replenishment  of  soil  moisture  and  recharge  of  groundwater,  and  eventually  in 
interflow  and  baseflow. 
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The  three  components  of  streamflow  have  been  treated  at  various  levels  of 
mathematical  sophistication  but  in  virtual  isolation  with  one  another.  Surface  flow 
has  been  studied  for  over  half  a  century  (Woolhiser,  1982),  and  as  a  result,  it  is 
understood  reasonably  well  (Hall,  1982).  Likewise,  baseflow  contribution  to 
streamflow  has  been  studied  for  nearly  30  years  and  it  too  is  understood  reasonably 
well  (Hall,  1982).  The  same,  of  course,  cannot  be  said  about  interflow.  This  is  not 
even  well  defined  and  is  least  understood.  Also,  least  understood  are  the  dynamic 
interactions  prevailing  amongst  these  components. 

Although  considerable  progress  has  been  made  in  mathematical  and 
numerical  treatment  of  the  boundary  value  problems  dealing  with  flows  over 
impervious  beds,  the  understanding  of  surface  flows  over  porous  beds  which 
dynamically  interact  with  subsurface  flow  is  quite  limited.  The  importance  of  this 
interaction  has  been  pointed  out  in  the  past  in  the  context  of  border  irrigation 
(Parlange,  1973),  and  in  the  study  of  flood  waves  in  ephemeral  streams  (Smith, 
1972).  These  studies,  however,  are  not  based  upon  a  coupled  set  of  equations 
pertaining  to  surface  and  subsurface  flow;  rather  the  attenuation  in  surface  flows  is 
included  by  considering  certain  infiltration  rate  with  time  lag.  The  dynamic  diffusion 
due  to  infiltration,  therefore,  remains  unaccounted  for. 

Freeze  (1972)  was  probably  the  first  to  develop  a  comprehensive  quantitative 
treatment  of  hillslope  hydrology  considering  explicit  interactions  between  near¬ 
surface  groundwater  flow,  surface  runoff  and  rainfall  intensity  patterns.  Rather 
limited  work  has  since  been  done  along  the  lines  of  Freeze.  Some  notable 
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examples  are  the  conceptual  model  of  Seven  and  Kirby  (1979),  the  model  of  Hillel 
and  Hornberger  (1979),  and  the  finite  element  model  of  Seven  (1977). 

The  most  recent  work  representing  a  major  step  forward  in  developing  an 
analytical  treatment  of  interdependent  surface  and  subsurface  hydrologic  processes 
is  by  Smith  and  Hebbert  (1983).  In  their  model,  the  hillslope  was  considered  to 
consist  of  two  soil  layers  with  the  lower  soil  capable  of  restricting  vertical  flow  at  the 
interface  creating  a  perched  aquifer  and  subsurface  stormflow.  Unsaturated  vertical 
flow  was  routed  by  a  kinematic  wave  method  and  linked  with  an  analytical  infiltration 
model.  Thus,  this  model  attempted  to  integrate  most  of  the  major  hydrologic 
response  mechanisms  presently  identified  as  contributing  to  the  hydrology  of  a 
simple  hillslope.  Other  hillslope  hydrological  models  (Gundy,  et  al.,  1985;  Stagnitti, 
et  al.,  1986)  and  surface  irrigation  models  (Walker  and  Humpherys,  1983;  Stagnitti, 
et  al.,  1986)  and  surface  irrigation  models  (Walker  and  Humpherys,  1983;  Ram,  et 
al.,  1983)  have  also  been  developed.  However,  none  of  these  models  developed  a 
method  to  compute  infiltration  rate  dynamically,  although  it  is  one  of  the  major 
factors  affecting  runoff  (or  advance  front)  and  surface  water  profile.  Most  of  the 
models  utilized  empirical  formulae  such  as  Kostiakov's  or  Green  and  Ampt's,  etc. 
Therefore,  the  prevailing  dynamic  process  between  surface  and  subsurface  flows 
remained  unaccounted  for. 

Toward  the  goal  of  eventually  accomplishing  a  continuum  model  for 
streamflow  synthesis,  three  related  areas  were  investigated:  (1)  subsurface 
unsaturated  flow,  (2)  flood  wave  propagation,  and  (3)  comparative  assessment  of 
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different  dynamic  wave  representations  of  shallow  water  wave  theory.  In  particular, 
we  focused  on  (a)  comparative  evaluation  of  the  kinematic-wave,  diffusion-wave, 
and  dynamic-wave  representations  of  the  shallow  water  wave  theory  ubiquitously 
applied  to  modeling  surface  runoff,  (b)  development  of  a  systems-  based  model  for 
infiltration,  and  (c)  modeling  movement  of  soil  moisture. 

2.  SURFACE  RUNOFF  MODELING 

The  surface  runoff  hydrology  was  investigated  along  three  lines:  (a) 
development  of  a  theory  of  errors  for  comparative  assessment  of  three  shallow 
water-wave  representations:  (i)  kinematic-wave,  (ii)  diffusion-wave,  and  (iii)  dynamic- 
wave;  (b)  development  of  discrete  linear  models  for  watershed  runoff;  and  (c) 
physically-based  Muskingum  methods  of  channel-flow  routing.  A  short  discussion  of 
each  is  in  order. 

2.1  Theory  of  Errors 

A  wide  range  of  problems  involving  free-surface  flows  can  be  modeled  using 
the  shallow  water-wave  (SWW)  theory.  The  SWW  theory  is  described  by  the  St. 
Venant  equations  or  their  approximations.  The  three  most  popular  representations 
of  the  SWW  theory  are  the  kinematic-wave  (KW),  diffusion-wave  (DW),  and  the 
dynamic-wave  (DYW)  approximations. 
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One  of  the  fundamental  questions  to  be  addressed  in  physically-based 
modeling  of  watershed  runoff  is  one  of  determining  the  appropriate  approximation  of 
the  shallow  water-wave  (SWW)  theory.  Of  the  different  approximations  of  the  SWW 
theory,  the  two  most  popular  approximations  are  the  kinematic-wave  (KW)  theory 
and  the  diffusion-wave  (DW)  theory.  How  accurate  are  these  approximations? 

Which  approximation  should  be  used  and  under  what  conditions?  What  is  the 
spatial  or  temporal  distribution  of  error  of  a  given  approximation?  What  is  the 
criterion  to  choose  between  these  approximations?  The  past  studies  have  dealt  with 
development  of  criteria  for  judging  the  adequacy  of  these  approximations.  However, 
these  criteria  are  point  values  and  do  not  relate  to  errors  resulting  from  use  of  these 
approximations.  Consequently,  the  error  in  space  and/or  time  is  not  known. 

The  larger  goal  of  this  study  was  to  develop  a  theory  of  errors  for  quantitative 
evaluation  of  the  adequacy  of  these  approximations,  and,  in  turn,  of  the  shallow- 
water-wave  theory.  However,  because  the  SWW  theory  consists  of  a  system  of 
nonlinear  partial  differential  equations  of  hyperbolic  type,  derivation  of  error 
equations  is  unattainable  at  this  stage.  Therefore,  some  realistic  simplifications 
were  made.  The  first  was  the  simplification  of  flows  being  time-independent. 

Steady  state  flows  are  ubiquitous  in  nature,  and  much  of  the  early  hydraulics  was 
based  on  this  assumption. 

Because  spatially  distributed  data  are  seldom  available,  it  was  assumed  that 
the  full  form  of  the  SWW  theory  or  the  dynamic-wave  representation  was  the  true 
representation  or  model,  and  was  capable  of  mimicking  the  behavior  of  the  real 
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world,  prototype  system,  and  the  kinematic-wave  and  diffusion-wave  approximations 
were  reasonable  approximations,  but  were  germane  to  conceptual  error.  The 
adequacy  of  these  approximations  is  well  documented  in  hydraulic  literature. 
However,  what  is  not  known  is  the  actual  error  and  its  distribution  in  time  and/or 
space,  as  well  as  its  relationship  to  flow  characteristics,  system  properties  and  initial 
and  boundary  conditions. 

The  theory  of  errors  can  serve  as  an  objective  criterion  for  judging  the 
adequacy  of  the  KW  and  DW  approximations  by  comparison  with  the  DYW 
approximation.  For  time-independent  flows,  the  theory  yields  error  as  a  function  of 
space  involving  infiltration,  and  boundary  conditions.  The  error  differential  equation 
is  ordinary  in  place  of  partial,  and  is  more  amenable  to  numerical  solution.  Even 
with  this  simplification,  analytical  solutions  are  not  possible  but  the  numerical 
solutions  are  much  simpler  and  easy  to  graph. 

Different  criteria  have  been  developed  to  evaluate  the  adequacy  of  the  KW 
and  DW  theories,  but  no  explicit  relations  either  in  time  or  in  space  between  these 
criteria  and  the  errors  resulting  from  these  approximations  have  yet  been  derived. 
Furthermore,  when  synthesizing  streamflow,  it  is  not  clear  if  the  kinematic-wave  and 
the  diffusion-wave  approximations  are  valid,  on  one  hand,  for  the  entire  hydrograph 
or  for  a  portion  thereof,  and  on  the  other  hand,  for  the  entire  channel  reach  or  for  a 
portion  thereof.  To  put  differently,  all  of  these  criteria  take  on  fixed  point-values  for 
a  given  rainfall-runoff  event.  This  study,  under  simplified  conditions,  derived  error 
equations  for  the  kinematic-wave  and  diffusion-wave  approximations  for  space- 
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independent  as  well  as  for  time-independent  flows.  These  equations  provided  a 
continuous  description  of  error  in  the  streamflow  hydrograph. 

With  these  considerations  in  mind,  errors  of  kinematic-wave  and  diffusion 
wave  approximations  were  derived  for  steady-state  channel  flows  subject  to  finite 
flow  at  the  upstream  end.  The  diffusion-wave  approximation  was  in  excellent 
agreement  with  the  dynamic  wave  representation  for  a  range  of  the  values  of  the 
Froude  number  and  the  kinematic-wave  number.  The  kinematic-wave  approximation 
was  also  in  good  agreement  with  the  dynamic  wave  representation,  but  for  a  limited 
range  of  the  values  of  the  Froude  number  and  the  kinematic-wave  number.  On  the 
other  hand,  the  approximate  diffusion-wave  analogy,  although  leading  to  analytical 
solutions,  was  not  accurate  and  should  not  be  employed. 

Under  two  different  initial  conditions  and  two  boundary  conditions,  solutions  of 
the  kinematic-wave  and  diffusion-wave  equations  were  derived  under  the 
simplification  that  the  flow  was  temporally  independent.  Thereafter,  error  equations 
for  the  KW  and  DW  theories  were  derived.  It  was  found  that  the  DW  theory  was 
quite  accurate  and  for  Froude  number  (Fo)  less  than  2  and  the  kinematic  wave 
number  (K)  greater  than  10,  the  DW  theory  would  be  an  accurate  representation  of 
the  SWW  theory.  Under  the  condition  where  there  was  no  downstream  control,  the 
KW  theory  was  an  accurate  representation  for  K  >  30,  K  Fq^  >  5  .  The  KW  theory 
does  not  accommodate  a  downstream  control  and  hence,  as  expected,  did  not 
accurately  represent  the  SWW  theory  for  the  entire  channel  under  any  of  the  two 
boundary  conditions.  Details  of  this  work  are  contained  in  Singh,  Aravamuthan  and 
Joseph  (1994). 
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For  space-independent  flows,  a  dimensionless  parameter  was  defined, 
reflecting  the  effect  of  initial  depth  of  flow,  channel-bed  slope,  lateral  inflow, 
acceleration  due  to  gravity,  and  channel  roughness.  For  time-independent  flows,  the 
dimensionless  parameter  was  the  product  of  the  kinematic-wave  number  and  the 
square  of  the  Froude  number.  By  comparing  the  kinematic-wave  and  diffusion  wave 
solutions  with  the  dynamic-wave  solution,  error  equations  were  derived  in  terms  of 
the  aforementioned  dimensionless  parameters.  The  error  equations  for  space- 
independent  flows  turned  out  to  have  the  form  of  the  Riccati  equation.  The  work  is 
described  in  Singh,  Aravamuthan  and  Joseph  (1993). 

2.2  Watershed  Runoff 

Discrete  linear  models  were  developed  for  estimating  runoff  and  sediment 
discharge  hydrographs  from  agricultural  watersheds.  A  regression  equation  was 
also  established  relating  runoff  rate  and  sediment  discharge.  Tested  on  five  small 
basins,  the  results  were  in  good  agreement  with  observations.  For  the  discrete 
linear  transfer  runoff  model,  the  values  of  the  integral  square  error  (ISE)  were 
generally  less  than  1%  for  all  calibration  events,  and  around  10%  with  the  average 
value  of  9.36%  for  all  verification  events.  For  the  discrete  linear  transfer  sediment 
model,  the  calibration  coefficient  of  determination  R  for  all  five  basins  was  more 
than  97%,  and  the  verification  R  was  more  than  91%  with  an  average  of  94.3%. 
Details  of  this  work  are  described  in  Wang  et  al.  (1991). 
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2.3.  Physically  -  Based  Muskingum  Method 

Flow  routing  in  channels  was  investigated  using  the  kinematic  wave  theory  as  well 
as  the  Muskingum  method.  For  the  latter  method  parameters  were  derived  from  the 
St.  Venant  equations.  Preliminary  testing  showed  that  this  method  of  parameter 
estimation  made  the  Muskingum  method  more  accurate  than  any  of  the 
conventional  methods.  The  kinematic  wave  method  was  investigated  for  perennial 
as  well  as  ephemeral  streams.  This  method  can  be  extended  to  include  flood  wave 
propagation  due  to  dam  rupture.  This  work  is  more  fully  described  in  Wang  and 
Singh  (1992). 

3.  SUBSURFACE  FLOW 

Modeling  of  flow  of  water  in  the  unsaturated  zone  is  far  from  complete, 
especially  at  the  field  scale.  Two  lines  of  inquiry  were,  therefore,  launched.  First,  a 
systerhs  approach  was  developed  to  model  infiltration  and  soil  moisture,  which  holds 
promise  for  unifying  different  infiltration  models  reported  in  the  literature.  This 
approach  can  also  relate  parameters  of  one  infiltration  to  another.  The  second  type 
of  approach  pertained  to  application  of  the  kinematic  wave  theory.  This  approach 
has  the  advantage  of  coupling  plant  root  extraction  and  evapotranspiration  with  soil 
water  dynamics. 

3.1  Infiltration  Modeling 

A  general  infiltration  model  was  derived  using  systems  approach.  The 
models  of  Horton,  Kostiakov,  Overton,  Green  and  Ampt,  and  Philip  are  some  of  the 
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example  models  which  are  shown  as  special  cases  of  the  general  model.  An 
equivalence  between  the  Green-Ampt  model  and  the  Philip  two-term  model  was 
shown.  The  general  model  also  provides  a  solution  for  the  Holtan  model  expressing 
infiltration  as  a  function  of  time.  This  solution  of  the  Holtan  model  does  not  appear 
to  have  been  reported  in  the  literature.  A  first-order  analysis  was  performed  to 
quantify  the  uncertainty  involved  with  the  generalized  model.  The  general  infiltration 
model  contains  five  parameters.  Two  of  the  parameters  are  physically  based  and 
can  therefore  be  estimated  from  the  knowledge  of  soil  properties,  antecedent  soil 
moisture  conditions,  and  infiltration  measurements.  The  remaining  three  parameters 
can  be  determined  using  the  least  squares  method.  The  model  was  verified  using 
ten  observed  infiltration  data  sets.  Agreement  between  observed  and  computed 
infiltration  was  quite  good.  This  work  is  more  fully  described  in  Singh  and  Yu 
(1990). 

3.2  Movement  Of  Soil  Moisture 

The  unsaturated  subsurface  flow  serves  as  a  link  between  surface  runoff  and 
groundwater  runoff,  and,  therefore,  occupies  the  central  position  in  the  streamflow 
continuum.  The  unsaturated  flow  region  is  the  upper  soil  matrix  which  also  is  the 
principal  source  of  interflow.  For  surface  flow,  infiltration  is  the  major  sink,  and  for 
groundwater  flow,  soil  moisture  percolation  is  the  principal  source  or  recharge. 

Much  of  the  mathematical  treatment  of  unsaturated  flow,  reported  to  date,  is  based 
on  the  Fokker-Planck  equation  or  Richards  equation.  Both  these  equations  are  of 
the  diffusion  type,  and  do  not  lend  themselves  to  analytical  solutions,  except  for 
overly  simplified  cases. 


If  one  ignores  the  effect  of  diffusion  and  models  soil  moisture  movement 
using  the  kinematics  wave  theory,  then,  under  certain  simplifying  but  realistic 
conditions,  it  is  possible  to  derive  analytical  solutions.  This  premise  was  pursued  in 
this  project.  Currently  available  one-dimensional  kinematic-wave  models  assume 
absence  of  sink  terms.  In  other  words,  once  the  water  gets  infiltrated,  it  is  either 
retained  by  the  soil  or  moves  downward  to  recharge  the  groundwater.  This 
assumption  is  not  tenable,  especially  in  agricultural  or  forest  watersheds.  In  this 
project,  an  effort  was  made  to  include  a  sink  term  in  modeling  of  soil  moisture.  This 
sink  term  may  represent  removal  of  soil  moisture  by  vegetation  through  the  process 
of  evapotraspiration. 

Recognizing  the  difficulties  of  modeling  unsaturated  flow  using  the  Richards 
equation,  a  novel  approach  was  developed  in  this  project.  This  approach  was 
based  on  the  kinematic  wave  theory  in  which  a  unique  relation  is  hypothesized 
between  the  flux  and  the  flow  concentration  or  the  hydraulic  head. 

The  fundamental  assumption  underlying  this  theory  is  that  the  moisture 
movement  is  gravity-dominated  and  the  hydraulic  conductivity-soil  moisture 
relationship  is  single-valued,  i.e.,  it  does  not  experience  any  hysteretic  effects. 
Although  these  assumptions  are  not  strictly  valid,  they  do  provide  a  reasonably 
accurate  approximation.  Another  advantage  of  the  theory  is  its  simplicity  and  that  it 
allows  analytical  solutions  under  simplified  conditions.  Under  more  complex 
conditions,  numerical  solutions  are  easily  derived. 
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This  unique  relation,  when  coupled  with  the  continuity  equation  expressing 
the  conservation  of  mass,  gives  rise  to  a  first  order,  nonlinear  hyperbolic  equation. 
Under  simplified  initial  and  boundary  conditions,  analytical  solutions  of  this 
kinematic-wave  equation  are  tractable.  Following  this  tract,  the  soil  moisture 
movement  was  modeled  with  the  use  of  the  kinematic-wave  theory.  However, 
attention  is  to  be  focused  on  certain  aspects  of  the  kinematic-wave  modeling  that 
are  not  apparent  at  the  first  glance.  First,  because  the  time  history  of  the  moisture 
front  distinguishing  between  wet  and  dry  soil  is  unknown,  the  kinematic-wave 
formulation  of  soil  moisture  movement  results  in  a  free-boundary  problem.  Second, 
natural  watersheds  have  vegetation  either  seasonally  or  throughout  the  year. 
Inclusion  of  vegetation  in  modeling  of  soil-moisture  movement  gives  rise  to  an 
additional  free  boundary,  further  complicating  the  model. 

With  these  considerations  in  mind,  a  kinematic-wave  model  was  developed 
for  simulating  the  movement  of  soil  moisture  in  unsaturated  soils  with  plants.  The 
model  involved  three  free  boundaries.  Analytical  solutions  were  derived  when  the 
plant-root  extraction  of  moisture  was  at  a  constant  rate,  and  the  upstream  boundary 
condition  was  time  independent.  If  these  assumptions  are  waived,  then  numerical 
solution  is  the  only  resort. 

With  the  use  of  this  theory,  a  comprehensive  analytical  treatment  has  been 
developed  for  soil  moisture  movement  with  plant-root  extraction.  The  treatment 
involves  free  boundaries  and  to  our  knowledge  this  has  not  been  dealt  with  in  the 
literature  so  far.  This  work  is  described  in  Singh  and  Joseph  (1994). 
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Abstract.  A  kinematic-wave  model  is  developed  for  simu¬ 
lating  the  movement  of  soil  moisture  in  unsaturated  soils 
with  plants.  The  model  involves  three  free  boundaries. 
Analytical  solutions  are  derived  when  the  plant  roots  are 
assumed  to  extract  moisture  at  a  constant  rate  and  the 
upstream  boundary  condition  is  independent  of  time. 
Numerical  solutions  are  the  only  resort  when  the  mois¬ 
ture  extraction  and  the  upstream  boundary  condition 
both  depend  on  lime. 


Much  of  the  mathematical  treatment  of  flow  in  unsatu¬ 
rated  porous  media  has  dealt  with  capillary-induced  flow 
(Smith  1983).  However,  there  exists  a  multitude  of  cases 
where  gravity  dominates  vertical  movement  of  soil  mois¬ 
ture.  Some  examples  include:  drainage  following  infiltra¬ 
tion,  the  water  percolation  deeper  into  the  soil,  the  verti¬ 
cal  movement  of  moisture  in  relatively  porous  soils  when 
rainfall  or  surface  fluxes  are  typically  of  the  order  of,  or 
less  than,  the  soil-saturated  hydraulic-conductivity,  to 
name  but  a  few.  For  treatment  in  such  cases,  the 
kinematic-wave  theory  is  simple  yet  reasonably  accurate. 

Although  Sisson  et  al.  (1980)  applied  the  kinematic- 
wave  theory  to  internal  drainage.  Smith  (1983)  was  prob¬ 
ably  the  first  to  apply  the  theory  to  develop  a  complete 
kinematic-wave  model  for  soil  moisture  movement. 
Charbeneau  (1984)  extended  Smith’s  work  to  solute 
transport,  and  Charbeneau  et  al.  (1989)  to  multiple  solute 
transport.  In  a  series  of  papers,  Germann  and  coworkers 
(Germann  1985;  Germann  and  Beven  1985,  1986;  Ger¬ 
man  et  al-  1987)  extended  the  application  of  the  theory  to 
infiltration  and  drainage  into  and  from  soil  macropores, 
as  well  as  to  microbial  transport.  Yamada  and  Kobayashi 
(1988)  discussed  the  kinematic  wave  characteristics  of 
vertical  infiltration  and  soil  moisture,  with  the  aid  of  field 
observations  on  tracers.  They  concluded  that  the  vertical 
infiltration  of  soil  moisture  had  the  characteristics  of 
kinematic  waves. 


Correspondence  to  V.  P  Singh 


In  this  Study,  a  kinematic-wave  model  is  formulated 
for  soil-moisture  movement  with  plant-root  extraction. 
Analytical  solutions  are  derived  under  the  condition  that 
the  moisture  extraction  by  plants  and  the  surface  flux  are 
both  constant  in  time.  This  condition  is  severe,  but  leads 
to  intriguing  and  useful  results.  For  more  realistic  con¬ 
ditions  numerical  solutions  are  the  only  resort,  these 
together  with  field  verification,  will  be  reported  in  the 
near  future. 

Kinematic-wave  model 

Formulation 

For  vertical  unsaturated  flow  with  plant-root  extraction 
and  with  the  neglect  of  the  flow  of  air,  the  governing 
equations  are  the  law  of  conservation  of  mass  and  a  flux 
law.  The  one-dimensional  conservation  of  mass  equation 
or  the  continuity  equation  can  be  expressed  as 


where  6  is  the  volumetric  moisture  content  (dimension¬ 
less,  volume  of  water  volume  of  soil),  q  is  the  vertical  flux 
of  soil  moisture  (L,T).  e(z.  x)  is  the  rate  of  soil  moisture 
loss  due  to  plant-root  extraction  at  time  r  (1  7).  r  is  the 
depth  (L)  measured  downward  from  the  ground  surface. 
T  is  the  time  (7)  that  water  has  stayed  at  location 
T  =  f  — vv(r),  vv(r)  is  the  time  history  of  the  moisture  ad¬ 
vance  front,  and  t  is  time  (7).  At  a  fixed  x,  eir)  is  essen¬ 
tially  the  rate  of  plant-root  extraction  of  soil  moisture  on 
a  unit  depth  basis  defined  as 


n 

where  n  is  the  soil  porosity,  and  r(f)  is  the  soil-moisture 
extraction  rate  per  unit  depth  at  lime  f  (1  7).  It  may. 
however,  be  reasonable  to  assume  ^(r.  r)  =  e(T).  This  as¬ 
sumption  implies  spatially  uniform  extraction,  that  may 
be  true  if  the  plant  roots  are  uniformly  developed  in  the 
soil. 
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The  flux  law  for  a  kinematic-wave  model  can  be  ex¬ 
pressed  as 

q{9)^K{e)  (3) 

where  K{0)  is  the  unsaturated  hydraulic  conductivity 
(L/T)  as  a  function  of  the  soil-moisture  content.  A  vari¬ 
ety  of  expressions  have  been  proposed  for  K  {9)  (Hsu  and 
Liu  1990).  Of  all  the  expressions,  the  best-known  perhaps 
is  the  Brooks-Corey  relation  (Brooks  and  Corey  1964) 
which  can  be  written  as 


I-* 

where  K,  is  the  saturated  hydraulic  conductivity  (LIT), 
9q  is  the  residual  value  of  9  below  which  water  cannot  be 
extracted  (or  moved)  by  capillary  forces,  9,  is  the  saturat¬ 
ed  water  content  =  porosity  =  n,  and  u  is  a  parameter 
typically  between  3  and  4,  and  is  related  to  the  pore-size 
distribution  index. 

In  Eq.(l),  T  is  unknown  and  a  relation  is  needed 
whereby  t  can  be  determined.  With  such  a  relation,  the 
formulation  of  the  kinematic-wave  model  will  be  com¬ 
plete.  To  that  end,  when  water  is  applied  to  the  soil  at  its 
surface,  it  infiltrates  into  the  soil  and  advances  vertically 
downward,  increasing  the  moisture  content.  If  the  soil  is 
initially  dry,  the  advance  front,  also  called  the  shock 
front,  will  define  the  interface  between  moist  and  dry 
parts  of  the  soil  matrix.  Let  r  =  s(r)  or  t  =  w(2)  be  the 
time  history  of  that  advancing  front;  this  time  history 
gives  the  space-time  coordinates  of  the  shock  front.  The 
front  is  a  free  boundary  and  has  to  be  determined  along 
with  the  solution.  An  equation  for  r  =  s(f)  or  t  =  w(r), 
therefore,  has  to  be  formulated.  This  can  be  accom¬ 
plished  by  observing  that 


u{9)  =  u(c,  r)  = 


qi9) 


KjO) 

9 


(5) 


where  u(0)  is  the  average  velocity  with  which  9  moves. 
Thus,  the  free-boundary  equation  is  obtained  by  replac¬ 
ing  r  by  sU)  or  r  by  w(r)  in  u(r,  t)  as 


^5(f) 

d! 

or 

dw{z) 

~dz~ 


=  u{s(r),  t)  = 


K{9{s{t),  r)) 
^(s(r),  t) 


s(0)  =  0 


=  [u(r,w(r))]‘^  = 


9{z,  Mz)) 
Kidiz,  w(z)))’ 


w(0)  = 


0. 


(6) 

(7) 


Equation  (6)  is  valid  for  9{z,  t)  >  0.  Equations  (6)  and  (7) 
require  the  advance  front  to  move  with  the  speed  immedi¬ 
ately  behind  the  front. 

The  kinematic-wave  model  formulation  consists  of  a 
partial  differential  Eq.  (U,  an  algebraic  Eq.  (3),  and  an 
ordinary  differential  Eq.  (6)  or  (7)  with  two  unknowns 
Oil.  f)  and  5(f)  or  uif).  In  order  to  derive  solutions  of  these 
equations,  the  following  can  be  assumed: 

0(r,0)=/(r),  r>0,  (8) 

9{Q.  t)  =  ^(f),  0<  f  <  7  (9) 

=  9o,  f>7 


where  T  is  the  duration  for  which  moisture  is  applied  at 
the  upstream  boundary.  Note  that  ^(0)  =/{0)  = 
Equations  (1)  and  (3)  can  be  combined  to  produce 


60  dK{6)  W 

ar  de  dz 


-e(z,  t). 


(10) 


For  notational  simplicity,  let 


M(0)  = 


dK(d) 

de 


dO/dt 

dO/dz 


dz 


dt 


Id  =  constant 


(11) 


which  is  referred  to  as  the  mobility  of  water  in  soil  by 
Irmay  (1956);  this  is  analogous  to  celerity,  commonly 
used  in  open  channel  flows. 


Solution 

Equation  (10)  is  satisfied  in  the  domain  bounded  by  the 
positive  t  axis  and  the  curve  r  =  w(z)  or  z  =  sit).  Solution 
of  Eqs.  (10)  and  (6)  can  be  derived  by  using  the  method  of 
characteristics.  Accordingly,  one  can  choose  cr  as  the 
parameter  on  the  f-axis  and  t  as  the  parameter  along 
the  characteristics.  The  characteristics  originate  from  the 
r-axis  on  the  segment  0  <  f  <  7.  The  solution  domain  can 
be  divided  into  two  subdomains  Dj  and  £>2,  with  the 
characteristics  z(t,  7)  serving  as  the  dividing  line  (called 
the  bounding  characteristic),  as  shown  in  Fig.  1.  The 
characteristic  curves  z(r,  <t),  9{t,<T)  passing  through  the 
points  (0,  (7,  g{a))  in  the  (z,  t,  6)  space,  satisfy 


(12) 

=M(t,a),  r(ff,  <t)  =  0. 

(13) 

t 


Fir.  1.  Characicnsiic  soluiion  domain  for  unsaturaicd  flov^  with 
plant-root  exiracuon 
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r(f,  (t)  is  the  position  at  time  t  of  the  soil-moisture  content 
6  which  was  at  z  =  0  and  at  f  =  a;  d{t,  <t)  is  the  soii-mois- 
ture  content  at  time  r.  The  plant-root  extraction  rate  de¬ 
pends  upon  the  nature  of  vegetation  and  availability  of 
soil  moisture.  If  e{x)  is  constant,  which  is  not  entirely 
unreasonable,  especially  when  the  soil-moisture  is  not 
limiting  plant  root  extraction,  then  ^(r)  =  e. 

Depending  upon  the  nature  of  the  function  g{t),  three 
cases  can  be  distinguished:  (1)  ^(r)  =  constant,  (2)  dg{t)/dt 
<0.  and  (3)  dg{t)/dt>0.  To  describe  the  variation  of  soil 
moisture  in  time  at  a  fixed  z  or  in  space  at  a  fixed  r,  one 
may  determine  dO/dz  and  d9/dt. 


Solution  for  domain 

This  domain  is  bounded  by  0  <  r  <  T,  z(f,  7),  and 
z  =  s(r).  Solution  of  Eq.(12)  is 

e(t,e)  g{<7)  -(t  -  (T)e  .  (14) 

Substitution  of  Eq.  (14)  in  Eq.  (13),  and  then  integration 
yields 

z(r,  a)  =  f  M(g{a)-{t  -a)e)dt .  (15) 

a 

With  the  use  of  Eq.  (4), 


M(0)=- 


K,a 


0-60^°’' 


{6,-60) 

Substitution  of  Eq.  (14)  in  Eq.  (16)  yields 

With  the  use  of  Eq.  (17),  Eq.  (15)  leads  to 


(16) 


(17) 


'g{a)  -  6o' 

a 

'g(a)-{t  -  ff)e  -6o 

e  1 

L  e.-6o  j 

t, 

1 

o 

(18) 


Equations  (14)  and  (18)  constitute  the  characteristic  solu¬ 
tion  for  domain  Dj .  By  eliminating  o  between  these  equa¬ 
tions,  0  can  be  expressed  as  a  function  of  z  and  t.  To  do 
that.  z(r.  c7)  must  be,  for  fixed  f,  cither  an  increasing  func¬ 
tion  of  (T  or  a  decreasing  function  of  o.  Differentiating 
Eq.  (18)  with  respect  to  a  yields 


z^(r,  a)  = 


5z(f,  a) 
da 


or 


cr)  = 


{O.-Oof 


\g'ia)[{g{a)-6Qf^^ -{Sit. 

^e{e{t.a)-eof]  .  (20) 


It  is  seen  from  Eq.  (20)  or  (19)  that  g  ia)  =  dg/da  <  0  is  a 
sufficient  condition  for  z^(t,  a)  <  0.  The  condition 


^'(<7)<0  includes  the  case  of  principal  interest  to  us, 
namely,  g{a)  =  constant. 


Case!,  ^(f)  =  constant:  If  ^(tr)  =  then  Eq.  (18)  be¬ 
comes 


e 


[f  -  ^0 1°  -a)  e-  Oq 

I  _  05  —  00  _  _  ^5  “  ^0 


.(21) 


By  coupling  Eqs.  (41)  and  (21),  6  can  be  expressed  as  a 
function  of  z  and  t  as 


0(z)-0o+(^,-^o) 


ez]^ 

.Js-Ooi  kJ 


(22  a) 


It  is  seen  that  in  domain  Dj ,  0  is  a  function  of  z  and  does 
not  depend  on  t.  The  soil  moisture  profile  is  nonuniform 
but  steady-state.  In  Eq.  (22  a),  there  are  two  terms  within 
braces.  The  second  term  is  due  to  water  uptake  and  the 
first  term  is  due  to  the  boundary  flux.  If  e  =  0,  then 
0(z)  =  0„.  The  nonuniformity  of  soil  moisture  is  a  direct 
consequence  of  the  plant-root  extraction.  If  e  were  zero, 
then  0  would  be  independent  of  z.  Thus,  the  alteration  of 
the  soil  moisture  profile  is  caused  by  plant  roots.  Of 
course,  the  degree  of  alteration  depends  on  the  value  of  e. 
On  a  short-time  scale,  the  value  of  e  is  relatively  small  in 
field.  Equation  (22a)  can  be  written  with  use  of  Eq.  (4)  as 


0(z)  =  6o  +  (0,-0o)(-^) 


(22  b) 


In  many  practical  cases,  ez/K(0J  is  much  less  than  1,  and 
attenuation  of  the  0-profile,  as  a  result,  will  be  small. 
When  Eq.  (22  a)  is  inserted  into  Eq.  (4),  the  soil-moisture 
flux,  with  the  aid  of  Eq.  (3),  is  obtained: 


K{e)  =  q(z)  =  K, 


(23) 


Similarly,  the  average  velocity  of  the  soil-moisture  move¬ 
ment  is  obtained  by  substituting  Eq.  (22)  in  Eq.  (5), 


Similarly,  substitution  of  Eq.  (22)  in  Eq.  (16)  yields  the 
mobility  of  water: 


(25) 


The  behavior  of  the  soil  moisture  in  this  case  can  be  fully 
described.  From  Eq.  (22)  at  a  fixed  z,  z  >  0. 


(26) 


and  at  any  r,  0<  /  <  7, 


ez 


ak. 


(27) 
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Equation  (27)  shows  that  at  any  depth  within  0  re¬ 
mains  constant  in  time  and  at  any  time  it  decreases  with 
increasing  z  from  the  ground  surface.  The  gradient  of  6 
increases  with  increasing  z.  This  corresponds  to  the  case 
of  a  wetting  event. 


Case  2  dg{t)jdt  <  0:  In  this  case,  60/0z  and  06/6r  are  ob¬ 
tained  from  Eqs.  (14)  and  (18)  by  noting  that 


00  ^  00/0CT  00  00/0(T 

0z  0z/0cr  ’  cr  0r/0(7 

To  that  end,  with  g'{a)  =  dg((T)/d(7, 
00 


(28) 


(29) 


Therefore,  by  dividing  Eq.  (29)  by  Eq.  (30), 
_ e{g'  +  e)(0,  -  0o) 


0r 


(31) 

(32) 


aK, 


fp-^o] 

a  - 1 

9  ^  {g  +  e) 

1  _  “  ^0. 

Lve.-v  j 

It  should  be  noted  that  in  this  characteristic  method  of 
solution,  r  was  chosen  as  the  parameter  along  character¬ 
istics.  Therefore,  in  Eq.  (18),  z  is  the  dependent  variable, 
and  a  is  the  independent  variable.  In  order  to  obtain 
00/0r,  one  can  choose  t  as  the  dependent  variable  and  z  as 
the  parameter  along  characteristics.  The  procedure,  how¬ 
ever,  remains  the  same  in  both  cases.  Employment  of  t  as 
the  dependent  vanable  results  in 


Cu  e  ^ 

=  gig  -  00)'- '  [{g  -  Bof  -  {B,- ■*’"  (32a) 

CCT 


and 

^  =  1  _i (3 _ 0  )«- > -doY-^ (B- eo)1" • 
ccT  e  K, 

(32  b) 

By  dividing  Eq.  (32a)  by  Eq.  (32  b), 


- .  e  9(9  -  Bof-  ‘  [{g  -  BoY  -i^iB-  00)*]'  ■ 

= - .(33) 

'''  e  -  g'ig  -  ej  -'[{9-0oY-^{B.-B^)T  “ 

The  variation  of  0  with  z  for  a  fixed  t  can  be  analyzed 
from  Eq.  (32).  If  e  <  g'ia)  then  00/0Z  >  0;  that  means  that 
0  increases  with  increasing  z.  On  the  other  hand,  if 
e^g'ia).  then  0  is  independent  of  z.  This  is  also  the  case 
when  z  >  K^/e {[g - 6q)  (6^  - Ojy.  If  e  >  g{a)  then  0  de¬ 
creases  with  z  but  can  also  have  a  complex  variation.  At 
a  fixed  z.  the  variation  of  <t  can  be  analyzed  with  the  aid 


of  Eq.  (33).  If  g'{a)  =  0,  then  0  is  independent  of  r.  If 
g'{<r)<0,  then  0  decreases  with  time.  If  ^  =  g'ia),  then  0  is 


a  decreasing  function.  At 

z  =  K,{{g-do)!(d-do)Yle, 

then 

bd  . 

—  =  0  and 
dt 

(34) 

dd  e(g'+e){d,-  d„Y 

dz  aK.g'ig-doY"'  ' 

(35) 

Equation  (34)  shows  that  0  becomes  independent  of  t 
whereas  Eq.  (35)  shows  that  as  a  function  of  z,  0  decreases 
with  z. 

Case  3,  dg(t)/dt  >  0:  This  can  be  analyzed  from  Eqs.  (32) 
to  (35).  However,  Eq.  (20)  reveals  that  when  g'{<T)>0, 
z^(t,  a)  >0  and  this  case  may  entail  shock  formation.  In 
that  case,  the  analytical  treatment  becomes  un wieldly  and 
is  not  discussed  in  this  paper. 


Determination  of  advance  front  (free  boundary 
of  domain  D^) 

In  order  to  determine  the  shock  front  r  =  w(z),  it  is  conve¬ 
nient  to  express  it  in  terms  of  the  parameter  a  introduced 
above.  If  one  considers  the  characteristic  curve  z(f,  cr)  in 
the  (z,  f)  plane,  then  it  will  intersect  the  free  boundary  at 
the  point  denoted  by  {i{(T),rj{<r))  as  shown  in  Fig.  1. 
Therefore,  z  =  t  =  rj{a)  is  the  parametric  representa¬ 
tion  of  the  free  boundary  in  terms  of  the  parameter  a.  For 
the  case  ^(f)  =  0„  =  constant,  and  e(T)  =  e  =  constant, 
determination  of  the  free  boundary  is  relatively  simple. 
From  Eq.  (22),  with  f  replaced  by  r/(a). 


z{rf(a\  cr)  =  <J(cT) 


(36) 

1 

From  Eqs.  (6),  (14)  and  (4),  and  noting  that  6  does  not  fall 
below  do. 


r/g.-goV  ( 


B,-(t}{(7)-ff}e-Bo 

B.-B„ 


dsjt) 

dt 


rt'iff) 


K{6(Ha),  rjiff))) 


K,  re, 

iB.-Bo)  . 


B.-B„ 


(37) 


Equations (36) and  (37) can  be  solved  for (aland  rj{a).  To 
that  end.  Eq.  (36)  can  be  differentiated  to  yield 


aK,  r 
"  id,  -  do)  L 


B„  -  inia)  -  a)e-do 

B.-dr, 


iifia)  -  1) . 
(38) 


On  eliminating  c'la)  between  Eqs.  (24)  and  (25). 


a  —  1 
which  gives 


ri(0)  =  0 


(39) 


(40) 
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Similarly,  eliminating  rj'ia)  between  Eqs.  (37)  and  (38) 
leads  to 

(0,-0o)  LU.-^oj 


e.-dp 
Os- do 


which  gives 


Equations  (40)  and  (42)  constitute  the  parametric  repre¬ 
sentation  of  the  free  boundary.  By  eliminating  <7  between 
Eqs.  (46)  and  (42),  and  noting  that  c  corresponds  to  z  and 
t]  to  ts  In  terms  of  q, 

■■?4[ftr-sr7il 

where  qo- q(9^).  Equation  (43)  or  (44)  expresses  the  ad¬ 
vance  front.  It  is  seen  from  Eq.  (43)  that  the  coordinates 
of  the  point  to  which  the  free  boundary  can  extend  is 
given  by 

(45) 

To  complete  the  solution  for  domain  ,  the  intersec¬ 
tion  of  the  bounding  characteristic  z(r,  T)  and  the  free 
boundary  r  =  s(f)  is  to  be  determined.  Let  that  point  be 
P  =  (r*,r*).  This  point  can  be  determined  explicitly  by 
equating  Eq.  (21)  with  <7  -  T,  r  =  t*  and  z  =.r*  to  Eq.  (43) 
with  r  =  r*  and  f  =  r*;  that  results  in 

r.  =  T  .  (46) 

a  “  I 

^  ^  _  r  _  __iZL_Tl 

e  liOs-Oo)  (fl-l)(^,-^o)J  J' 


For  the  case  when  ^(f)  is  not  constant,  we  use  the  same 
parametric  representation  as  before.  Therefore,  analogous 
to  Eq.  (36), 

e{U^  —  aQ} 

(48) 

and  analogous  to  Eq.  (37), 


{9,  -  9J 


[g{(X)-[r]{c)-a]e-9Q] 


Differentiating  Eq.  (48)  and  then  substituting  the  dilTeren- 
lial  in  Eq.  (49)  yield 

^  __a _ g{a)[g{c)  -  ^or~* 

4cr'^e(fl-l)  [g{c)-{n{(T)-<T)e-9QY^^ 


Equation  (50)  cannot,  in  general,  be  solved  in  terms  of 
simple  forms  of  g(t),  an  analytical  solution  for  ?/(cr)  is 
tractable.  With  substitution  of  this  solution  in  Eq.  (48), 
q((T)  is  obtained.  If,  however,  g'((T)  =  0,  Eq.  (50)  reduces  to 
Eq.(39). 


Another  method  to  determine  the  advance  front 
The  velocity  of  the  front  can  be  expressed  as 


Representing  by  and  0,  and  taking  advantage 
of  Eqs.  (22)  and  (23),  Eq.  (46)  can  be  written  for  the  case 
g{t)  =  constant,  as 


ds{t) 

dt  ^9,-9o 


_  „  U/fl 


iq  -  ez)' 


which,  when  integrated,  leads  to 

'  e  e  [[kJ  a{6,-eo) 
which  is  the  same  as  Eq.  (44). 


Solution  for  domain  Dj 

Domain  D2  is  bounded  by  7  <  r,  r  =  r(f,  T),  z  =  s(f)  be¬ 
yond  the  point  p  =  (r«,  r*).  In  this  case,  ^(0,  r)  =  ^0  for 
t>  7;  9q<9^s  Let  g{t,  £)  be  a  continuous  function  of  t 
coinciding  with  ^(r)  on  0<  r  <  7,  ^(r,  £)  =  ^0  on  r  >  7  + 
as  shown  in  Fig.  2.  Then,  it  is  seen  from  Eqs.  ( 1 4)  and  (18) 
that  the  characteristics  r  =  r(r,  cr,  £),  7  <  a  <  T-hc,  cover 
the  region  above  2  =  2(r,  7)  completely. 


K,  f  9{t,  a,  £)  -f-  (f  - 

2(f,  ff,  £)  =  ^  ^  - - - — 


9U.  a,  e)  —  9q 


e{a  —  1 


■(g'{a)  +  e).  r;(0)  =  0. 


t  F'S2.  Function  (/(;.  t:) 
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By  letting  ^  0  and  a  -►  7,  one  obtains 


•  (54  b) 


This  equation  expresses  the  soil-moisture  0  as  a  function 
of  :  and  r  in  domain  9  varies  from  9^  to  9q.  The 
reduction  in  9  is  caused  by  the  plant-root  extraction,  and 
hence  the  water  uptake  is  the  dominant  factor  for  varia¬ 
tion  of  water  content.  This  can  be  easily  shown  for  a 
special  case  when  a  is  an  integer.  To  that  end,  Eq.  (54  b) 
can  be  written  as 


[9-9o 

\9s-9o 


\re{t^T) 

IL  ^-^0 


(55  a) 


With  use  of  the  bionomial  theorem,  this  can  be  written  as 


^-^0  Y 
o.-eoJ 


mr) 


(55  b) 


=  K, 


e-dp 

e.-d. 


a  (a-\ 

I 


{er 


(t-T) 

0-00 


It  is  seen  that  each  term  in  the  series  is  affected  by  the 
plant-root  extraction  e  raised  to  the  power  (a  — /  —  I). 
Assume  that  only  the  first  term  of  the  series  is  most  signif¬ 
icant.  Then. 


=(0.  r)  =  K, 


f0-0o 

Us -00 


a  e 


a  —  1 


r-  T 
9-9, 


(55c) 


Thu.«  the  reduction  in  9  is  primarily  caused  by  water 
uptake. 

The  characteristics  in  domain  Dj  are  given  by 


r  =  ^ 

1 

o 

1 

/e.-e{t-T)-0oy'\ 

e- 

1 

o 

1 

[  0S-0O  )} 

(56) 


where  is  the  initial  soil-moisture  9o<9i<9^,  that  is 
carried  out  by  the  i-lh  characteristic. 


Determination  of  free  boundaries  for  domain  Dj 


Domain  D,  has  two  free  boundaries  FB2  and  FB3  as 
shown  in  Fig.  3.  Along  FB3  that  starts  at  f  =  7  and  x  —  0, 
0{z,  t)  =  Of^-Ahus.  this  is  the  locus  defined  by  ^o-T'his  is  the 
time  history  of  the  drying  front  as  it  advances  from  x  =  0. 
From  Eq.  114).  we  have 

0q=  gicr)  —  it  —  o)  e.  7<(7<7  +  £.  (57) 

Thus,  from  Eq.  (18).  we  obtain 


r(f.  G,c)  =  — 
e 


(f  -  <7|  cl® 

9.  -  9o  ^ 


(58) 


Letting  £  0,  a  7.  we  get  the  free  boundary  FB3 


c(f-  7)1® 
9.  -9,  ^  • 


Between  r  =  0  and  Eq.  (59).  f^(r.  t)  —  0,. 


(59) 


t 


Fig.  3.  Characteristic  solution  domain  for  the  case  ^'(0  =  0, 
0^/<  7 


t 


Fig.  4.  Characteristic  solution  domain  for  the  case  when 
e  =  constant  >0.  —e<g'{!)<0.0<f<T 


It  IS  seen  from  Eq.  (56)  that  the  maximum  distance  that 
the  characteristics  will  travel  is  given  by 


which  is  the  same  as  given  by  Eq.  (45).  This  shows  that 
FB2  (the  locus  9{z,  f)  =  ri  a  vertical  wall  with  lower 
end  given  by  Eq.  (45)  and  the  upper  end  by  Eq.  (56)  and 


t*  =  T+-(0,-Ho).  (61) 

e 

In  the  case  when  t/(f)  is  not  constant,  we  gel  the  paramet¬ 
ric  representation  of  FB2; 


(62) 


MOISTURE  FLUX  MOISTURE  FLUX 
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0  10  20  30 


Fig.  7.  Moisture  flux  as  a 
function  of  time  for  fixed 
values  of  depth 


TIME.  H 


Fig.  8,  Moisture  flux  as  a 
function  of  depth  for  fixed 
values  of  time 


Total  moisture  content 

The  total  moisture  content  W  varies  within  the  region, 
and  can  be  determined  by  integrating  the  moisture  profile 
from  r  =  0  to  a  particular  depth  z  =  L.  Let  L  <  z*.  Refer¬ 
ring  to  Fig.  1.  at  f  =  0,  the  initial  total  moisture  content 
m;  =  OqL.  For  0  <  f  <  r.  say,  r  =  r, ,  when  g{t)  =  con¬ 
stant. 


where  z^  is  the  point  of  intersection  of  r  =  r  j  and  r  =  s(f). 
For  t  >  f*. 

W{t.  L)  =  I  t)dz  (66) 

0 

where  f)  is  prescribed  by  Eq.  (55).  An  explicit  solution 
of  this  integral  is,  however,  not  tractable. 


H  =  r ,  -F  (if) (L  — ’  z , ) 


(65) 
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DRAINAGE  PERIOD  »  6  H.  e  *  0.2/DAY 
BROOK  AND  COREY  RELATION 
FIXED  DEPTH  *  10,  20.  30  CM 


Fig.  9.  Average  velocity  as  a 
function  of  time  for  fixed 
values  of  depth 


DRAINAGE  PERIOD  «  6  H.  6  «  0.2/OAY 
BROOK  AND  COREY  RELATION 
FIXED  TIME.  2,  7.  10H 


Fig.  10.  Average  velocity  as 
a  function  of  depth  for  fixed 
values  of  lime 


depth.  CM 


An  example  of  application 

A  wetting  event  with  moisture  applied  at  the  ground  sur¬ 
face  was  considered,  with  K{9)  given  by  Eq.  (4).  The  soil 
was  Glendale  clay  loam,  investigated  by  Sisson  et  al. 
(1980),  with  the  following  properties:  0,  =  0.52,  = 

0.246,  Kj  =  1  m/day,  and  a  =  4.25.  The  moisture  was 
applied  at  the  rate  of  2  cm/h  for  6  hours.  These  values 
were  used  by  Charbeneau  (1984).  The  plant-root  extrac¬ 
tion  e(t)  was  considered  constant  but  values  of  e  =  0.1, 
0.2,  0.3,  0.4.  and  0.5 /day  were  considered.  It  should  be 
emphasized  that  these  values  were  used  to  illustrate  the 
procedure  presented  in  the  paper.  If  n  =  0.52,  then 
c  =  0.1/day  would  correspond  to  r(r)  =  0.052/day.  For 
one  meter  of  soil  column,  this  would  translate  into 
5.2  mm  of  plant-root  extraction  for  a  day.  Likewise, 
e  =  0.5/day  would  yield  r{t)  =  0.25/day  or  25  mm/day. 
This  rate  is  clearly  high  but  the  objective  here  was  to 
intentionally  show  the  modulation  of  soil  moisture  profile 
by  the  plant-root  extraction. 

The  kinematic  wave  model  was  applied  to  compute 
0(r,  f),  q{z,t\  and  u(z,  r).  Then  ^(f),  u{t)  and  ^(f)  were 
obtained  at  selected  values  of  z.  Likewise,  d{z)  and  q{z) 
were  computed  at  selected  values  of  t.  For  a  sample  case 
of  e  =  0.2 /day,  computations  are  illustrated.  Figure  5 
shows  ^(r)  at  2  =  10,  20,  and  30  cm,  and  Fig.  6  shows  6(z) 
at  f  =  2,  7  and  10  h.  The  moisture  profile  has  realistic 
appearance.  The  moisture  flux  is  plotted  as  a  function  of 
f  in  Fig.  7,  and  as  a  function  of  z  in  Fig.  8.  The  average 
velocity  is  plotted  as  a  function  of  t  in  Fig.  9,  and  as 
function  of  r  in  Fig.  10.  Both  q  and  u  appear  to  be  realistic 
in  their  variation.  The  effect  of  increasing  e  was  to  com¬ 
press  domain  .  In  other  words,  after  the  moisture  flux 
was  ended  at  the  ground  surface,  the  moisture  profile 
returned  to  its  initial  state  more  quickly  where  the  mois¬ 
ture  extraction  rate  was  higher.  For  higher  e,  the  charac¬ 
teristics  in  domain  were  more  curved. 


Conclusions 

The  kinematic  wave  model  for  soil  moisture  movement 
has  been  formulated,  and  its  solution  domain  involves 
three  free  boundaries.  Analytical  solutions  are  tractable 
for  the  case  when  the  upstream  boundary  condition  is 
independent  of  time.  For  a  time-dependent  boundary 


condition,  numerical  solutions  are  the  only  resort.  In  this 
model,  the  advancing  front  of  water  or  wetting  front  is 
necessarily  an  advancing  shock  wave.  It  is  plausible  to 
assume  that  the  moisture  content  at  the  wetting  front  is 
Qq.  Since  the  kinematic  wave  theory  cannot  account  for 
this  assumption,  the  diffusion  wave  model  has  to  be  em¬ 
ployed  in  that  case. 
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for  time-independent  flows  in  infiltrating  channels 
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Abstract  Time-independent  (or  steady-state)  cases  of 
channel  flow  were  treated  and  errors  of  the  kinematic-wave 
and  diffusion-wave  approximations  derived  for  finite  flow 
at  the  upstream  end.  The  diffusion-wave  approximation 
was  found  to  be  in  excellent  agreement  with  the  dynamic 
wave  representation,  with  error  magnitudes  of  0.2%  for 
values  of  KFq>1  .5^  where  K  is  the  kinematic-wave  num¬ 
ber  and  Fq  is  the  Froude  number.  Even  for  small  values  of 
KFq  (e.g.,  KFq-0.15),  the  errors  were  typically  in  the 
range  of  1 .3  to  3.7%.  The  approximate  analytical  diffusion- 
wave  solution  performed  poorly  with  error  magnitudes 
greater  than  30%  even  for  large  values  of  KFq  .  The  kine- 
matic-w^ave  approximation  was  also  found  to  be  in  good 
agreement  with  the  dynamic-wave  representation  with  er¬ 
rors  of  about  1.2%  for  KFq=1,5  and  varying  from  15  to 
44%for  A:Fo^=0.75. 


In  overland  flows,  the  steady  state  is  attained  at  the  outlet 
for  constant  rainfall  after  the  depth  of  flow  reaches  the  equi¬ 
librium.  The  same  is  true  for  flow  in  a  channel  subject  to 
constant  lateral  inflow.  For  a  channel  receiving  a  constant 
inflow  of  long  duration  at  the  upstream  boundary,  the  flow 
at  the  downstream  end  would  reach  the  equilibrium  and 
this  frequently  occurs  in  border  and  furrow  irrigation.  De¬ 
spite  occurrences  of  steady  state  or  time-independent 
flows,  they  have  not  received  much  attention  in  hydrology 
(Morris  1978).  The  steady-state  solution  aids  in  under¬ 
standing  the  nature  of  the  water-surface  profile.  It  may  help 
determine  the  condition  for  use  of  zero  depth  in  place  of 
zero  influx  at  the  upstream  boundary.  When  the  rainfall  du- 
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ration  is  much  greater  than  the  time  of  equilibrium,  the 
steady-state  water  surface  profiles  are  very  useful. 

Mathematical  treatment  of  steady-state  flows  is  based 
on  the  shallow  water-wave  (SWW)  theory  or  its  simplifi¬ 
cations,  including  the  fcinematic-wave  (KW)  or  the  Jiffu- 
sion-wave  (DW)  approximation,  in  conjunction  with  ap¬ 
propriate  boundary  conditions.  A  comprehensive  discus¬ 
sion  of  steady-state  flows  using  the  <iiffusion-wave  (DW) 
approximation  was  presented  by  Govindaraju  et  al.  (1988) 
including  both  numerical  and  analytical  results  for  flux- 
type  boundary  conditions.  The  upstream  boundary  condi¬ 
tion  was  one  of  zero  inflow.  Both  the  zero  depth-gradient 
and  the  critical-depth  downstream  boundary  conditions 
were  investigated.  For  steep  slopes,  the  upstream  bound¬ 
ary  condition  of  zero  depth  was  found  to  be  justifiable.  It 
was  shown  that  the  critical-depth  condition  at  the  down¬ 
stream  boundary  was  a  stringent  condition  which  might 
give  rise  to  problems  for  certain  ranges  of  parameter  val¬ 
ues  when  seeking  a  numerical  solution.  Furthermore,  Go¬ 
vindaraju  et  al.  (1988)  employed  an  analytical  approxima¬ 
tion  in  the  form  of  a  cubic  approximation  of  the  DW  model 
for  the  zero  depth-gradient  downstream  boundary  condi¬ 
tion.  The  cubic  approximation  was  found  to  be  accurate 
when  Fq  K  was  sufficiently  small,  with  Fq  being  the  Froude 
number  and  K  the  kinematic  wave  number. 

Parlange  et  al.  (1989)  improved  the  approximation  of 
Govindaraju  et  al.  (1988)  by  using  a  Taylor  series  approx¬ 
imation,  and  showed  that  the  improved  approximation  and 
earlier  approximation  would  provide  the  upper  and 
lower  bounds  of  the  numerical  solution.  Most  studies  deal¬ 
ing  with  evaluating  the  adequacy  of  the  kinematic-wave 
(KW)  or  diffusion-wave  (DW)  approximation  have  dealt 
with  space-time  dependent  flows.  Criteria  forjudging  the 
adequacy  of  these  approximations  have  been  derived  as 
point  values  either  in  terms  of  K  or  KFq  for  given  flov'  sit¬ 
uations,  but  they  do  not  relate  to  errors  in  time  and/or  space. 
Moreover,  an  explicit  treatment  of  steady-state  flows  has 
not  been  included  in  these  studies. 

Pearson  (1989)  examined  the  criteria  for  using  the  KW 
approximation  to  the  St.  Venant  (SV)  equations  for  shal¬ 
low  water  flow.  For  steady-state  one-dimensional  flow 
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over  a  plane  he  showed  that  for  the  condition  of  zero  flow 
at  the  upstream  boundary,  the  assumed  steady-state  up¬ 
stream  depth  is  generally  nonzero,  but  is  implicitly  as¬ 
sumed  to  be  zero  in  the  kinematic-wave  approximation.  By 
plotting  contours  of  dimensionless  steady-state  upstream 
depth  over  a  2-parameter  plane  (Fq,  KF^)  he  obtained  a 
new  criterion  as  K>3-^5IFq  for  kinematic-wave  model¬ 
ing. 

Parlange  et  al.  (1990)  were  probably  the  first  to  under¬ 
take  an  investigation  of  errors  in  the  KW  and  DW  approx¬ 
imations  by  comparing  their  predictions  with  the  numeri¬ 
cal  solution  of  the  5t.  Venant  (SV)  equations  under  steady 
state  conditions.  It  was  shown  that  the  two  approximations 
could  have  significant  errors  even  for  critical  flow  and 
fairly  large  values  of  the  kinematic  wave  number.  When 
the  KW  approximation  was  inaccurate,  the  improvement 
by  the  DW  approximation  was  only  modest.  Parlange  et 
al.  (1990)  then  proposed  a  more  accurate  approximation 
when  the  kinematic  wave  number  was  large.  They  sug¬ 
gested  splitting  the  solution  of  the  SV  equations  in  two  re¬ 
gions,  one  near  the  downstream  end  of  the  plane  and  the 
other  covering  most  of  the  plane. 

However,  no  explicit  relations  as  a  function  of  space 
between  these  criteria  and  the  errors  resulting  from  the  KW 
or  DW  approximation  have  been  derived  as  yet.  Also,  these 
criteria  did  not  include  the  effect  of  infiltration  although 
this  is  vital  for  surface  irrigation.  As  a  result,  the  actual 
error  of  these  approximations  as  a  function  of  space  is 
usually  not  known.  Furthermore,  it  is  not  evident  in  hydro- 
logic  modeling  if  the  KW  and  DW  approximations  are  valid 
for  the  entire  length  of  the  channel  or  a  portion  thereof,  or 
for  infiltrating  channels. 

The  objective  of  this  study  was  to  derive,  under  simpli¬ 
fied  conditions,  these  errors  as  a  function  of  space  for  the 
KW  and  DW  approximations  for  time-independent  flows 
in  infiltrating  channels. 


Shallow  wateMvave  (SWW)  theory 
for  time-independent  flows 

Governing  equations 

The  SWW  theory  for  time-independent  flows  can  be  de¬ 
scribed  by  the  St.  Venant  (S  V)  equations.  For  a  time-inde¬ 
pendent  flow  in  an  infiltrating  channel  of  length  L,  these 
equations  can  be  written  in  one-dimensional  form  on  a  unit 
width  basis  as: 

Continuity  equation, 

Ao, /,)  =  -/  ,1, 

and  momentum  equation, 

=  (2) 

where  h  is  depth  of  flow  (L),  u  is  local  mean  velocity 
(L/r),/is  uniform  infiltration  rate  (L/T),  g  is  accelera¬ 


tion  due  to  gravity  {L/T^),  x  is  space  coordinate  in  the  di¬ 
rection  of  flow  (L),  Sq  is  bed  slope,  and  S/is  frictional  slope. 
Note  Q-u  h  is  discharge  {L?I{TL))  per  unit  width,  and  Eq. 
(2)  neglects  the  effect  of  momentum  exchange  between 
longitudinal  flow  and  infiltration.  This  amounts  to  assum¬ 
ing  that  infiltration  occurs  in  the  vertical  direction  only.  Sf 
can  be  approximated  as 

where  (T^/L)  is  some  resistance  parameter.  If  the  Chezy 
relation  is  used  then  where  C  is  Chezy ’s  resis¬ 

tance  parameter. 

Equations  (1)  and  (2)  are  the  governing  equations  for 
the  dynamic  >vave  (DYW)  representation  for  time-indepen¬ 
dent  flows.  The  KW  approximation  is  based  on  Eq.  ( 1 )  and 
Eq.  (2)  with  the  left  side  deleted, 

g{SQ~Sf)  or  (4) 

which  can  be  expressed  as  Eq.  (3).  The  DW  approxima¬ 
tion  uses  Eq.  (1)  and  Eq.  (2)  with  the  convective  acceler¬ 
ation  term  deleted. 


It  is  useful  to  express  the  SV  equations  in  dimension¬ 
less  form.  To  that  end,  the  following  normalizing  quan¬ 
tities  are  defined:  (2o  =  discharge  per  unit  width  (L^/T)  of 
the  channel  at  the  upstream  boundary;  //o  =  normal  depth 
(L)  of  flow  corresponding  to  Qq  at  the  upstream  boundary; 
C/o  =  normal  velocity  (L/T)  for  Qo(f2o/^o)L-^o= normal¬ 
izing  distance  (L)  computed  as  Qo/f.  where/  is  the  aver¬ 
age  rate  of /infiltration  during  the  period  of  irrigation 
obtained  as/=jQ/(w)  dw/T, /(f)  =  infiltration  rate  (L/T), 
T -  duration  of  irrigation,  Tq  =  normalizing  time  (T)  com¬ 
puted  as  To  =  //o//,  and  Fq  =  normalizing  Froude  number 
=  should  be  noted  that  the  quantities  quali¬ 

fied  as  “normal”  correspond  to  the  normal  flow,  whereas 
those  qualified  as  “normalizing”  are  employed  to  nondi- 
mensionalize  the  flow  variables.  With  the  use  of  the  above 
normalizing  quantities,  the  following  dimensionless  quan¬ 
tities  can  be  defined: 


Xq  Mq  Uo 

(6) 

o  -  Q  f  F  M* 

Qo  r  * 

(7) 

Recall  that  the  Froude  number  F  is 

F  =  Froude  number  =  — 

(8) 

and  the  kinematic  wave  number  K  is 

II 

(9) 

With  the  introduction  of  the  above  normalizing  quan¬ 
tities,  Eq.  (1)  becomes 
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djc. 


K  K)  =  -f* 


and  Eq.  (2)  becomes 
d«.  f?  dt. 


r 


«* '' 

^--r 
K ) 


(10) 


(11) 


Equations  (10)  and  (11)  are  the  governing  equations  for 
the  DYW  representation  for  time-independent  flows.  The 
KW  approximation  is  based  on  Eq.  (10)  and  Eq.  (11)  with 
the  left  side  deleted, 


K 


M*  ^ 

•"TT 

K) 


=  0  or  ui  =h^. 


(12) 


This  is  actually  the  formal  limit  of  Eq.  (11)  when  1/KFq 
0.  The  DW  approximation  uses  Eq.  (10)  and  Eq.  (11) 
with  the  convective-acceleration  term  deleted, 


Jl^ 
Fi  dr* 


=  K 


f- 


(13) 


This  can  be  deduced  from  Eq.  (1 1)  by  letting  Fq  ->  0  for 
fixed /sTFo^ 


Boundary  conditions 

The  boundary  condition  can  be  specified  at  the  upstream 
boundary  as  follows: 

fH0)=ho,  M(0)  =  Mo,  Q(0)=Qo.  (14) 

In  dimensionless  terms, 

h*iO)=ho^,  M*(0)=Mo*,  e*(0)  =  eo*-  (15) 


Definition  of  error 


The  relative  error  E  can  be  defined  by 


(16) 


where  is  the  dimensionless  solution  from  KW  or  DW 
approximation  and  Sp  is  the  dimensionless  solution  from 
DYW  representation.  The  solution  can  be  either  in  terms 
of  depth  ()i* ),  velocity  («* ),  or  discharge  (g* ).  Thus, 

ud  hp  Qp 

where  the  subscripts  K  and  D  correspond  to  the  KW  (or 
DW)  and  DYW  solutions,  respectively. 

The  differential  equation  for  error  can  be  obtained  by 
differentiating  Eq.  (16)  with  respect  to  as 

d£^(F  +  l)d5;,  (F-H)^  d5o 

dx*  5;,  dr*  5^  dr*  ' 


Note  that 


Sd  = 


E  +  \  ■ 


(19) 


Equation  ( 1 8)  was  used  to  define  the  error.  The  differen¬ 
tial  equation  for  error  can,  however,  be  defined  without  ex¬ 
plicitly  knowing  Sp  so  long  as  is  known,  as  seen  from 
Eq.  (19). 

For  the  above-specified  boundary  condition,  twelve 
cases  involving  Fo  =  0-1,  0.5,  1.0;  and  K=3,5,  10  and  30, 
were  considered  for  computing  errors.  The  KW,  DW  and 
DYW  solutions  were  computed  in  terms  of  depth,  veloc¬ 
ity,  and  discharge  for  all  of  these  cases. 


Hydrodynamic  solutions 

Kinematic-wave  solution 

Equation  (10),  subject  to  the  upstream  boundary  condition 
given  by  Eq.  (15),  has  the  solution: 

W*  ^*  =  2*=-/*^*  +  ^  ,  Uo^  =  Qo^=l  (20) 

where  a  is  the  dimensionless  discharge  at  the  upstream 
boundary  which  equals  one.  With  introduction  of  Eq.  (12), 

=  )"-'  =  (! (21) 

which  is  the  KW  solution.  When  expressed  in  terms  of  di¬ 
mensionless  depth,  it  becomes 

K  =  =  -/*  j:*  .  ( 22) 

Equations  (21)  and  (22)  show  that  the  flow  will  advance 
as  far  as  X^  =  a/f^  =  I//"*,  and  clearly  depends  on  the  up¬ 
stream  condition  and  the  rate  of  infiltration. 

The  KW  solution  (depth,  velocity,  and  discharge)  was 
computed  for  all  12  different  cases  involving  values  of  the 
Froude  number  (Fo=0.1,  0.5,  and  1.0),  kinematic  wave 
number  (K=3,  5,  10,  and  30),  lateral  inflow  ^*=0,  infil¬ 
tration  rate /* =0.5,  and  dimensionless  upstream  discharge 
(a=  1.0).  The  solution  in  terms  of  flow  depth  for  a  sample 
case  (Fo=0.5)  is  shown  in  Fig.  1.  The  KW  solution  is  not 
dependent  on  K  or  Fq.  The  depth  decreases  with  the  lon¬ 
gitudinal  distance.  This  is  due  to  the  fact  that  the  total  in¬ 
filtration  increases  from  upstream  to  downstream. 


Diffusion-wave  solution 


Equation  (10),  subject  to  the  condition  in  Eq.  (15),  has  the 
solution  given  by  Eq.  (20).  Substitution  for  from  Eq. 
(20)  into  Eq.  (13)  yields 


«-/*  -t* 


=  K  Fit 


1-- 


(«-/* 


(23) 


Differentiation  and  simplification  of  Eq.  (23)  produce  the 
DW  equation  in  terms  of  «* . 
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Kinematic  oia^^e  approximation  :  Froude  number 

Upstream  boundary  condition  :  Constant  discharge  =  1 
Lateral  in-flouj  q  =  0.00  Infiltration  F  =  0.50 

©«  K  =  3.00  K  =  5.00  -M-  K  =  10.00  >e<  K 


0.50 

30.00 


Fig.  1  Dimensionless  flow  depth  by  the  kinematic-wave  approxi¬ 
mation  as  a  function  of  dimensionless  distance  when  the  upstream 
discharge  is  1.0,  Froude  number=0.5.  and  is  variable 


^  ^  A'Fq-  ul-K  Fq-  ul  (1  (1  -/^  xj 


(24) 


Its  analytical  solution  is  not  tractable. 

The  DW  approximation  can  also  be  expressed  in  terms 
of  /j*  as 


1- 


hi 


(25) 


Equation  (25)  shows  that  the  depth  gradient  becomes  zero 
when 


A*  =  (l-/*-^*)^^  (26) 

which  is  the  KW  solution  given  by  Eq.  (22).  When 
/*  x*  =  1,  the  depth  gradient  equals 


(27) 


An  attempt  was  made  to  integrate  Eq.  (25)  by  prescrib¬ 
ing  the  boundary  condition  at  the  upstream  end,  but  this 
resulted  in  a  solution  which  was  not  realistic.  This  can  be 
attributed  to  the  highly  nonlinear  nature  of  Eq.  (25)  and  the 
presence  of  subcritical  flow.  A  viable  solution  can  be 


obtained  by  prescribing  a  downstream  boundary  condi¬ 
tion  and  marching  the  solution  from  downstream  to  up¬ 
stream. 

A  solution  could  not  be  obtained  for  the  downstream 
depth  h^=0.0.  A  very  small  but  finite  depth  A,.  =  0.04  was 
prescribed  at  the  downstream  end.  By  employing  this  pro¬ 
cedure,  the  diffusion  wave  solution  (dimensionless  head, 
velocity  and  discharge)  was  computed  for  all  12  cases 
involving  various  values  of  Fq  and  K.  For  a  sample  case 
(Fo=0.5),  the  solution  in  terms  of  dimensionless  depth  is 
shown  in  Fig.  2.  As  can  be  seen  from  the  figure  for  a  con¬ 
stant  downstream  depth  of  /i*=0.04,  the  upstream  depth 
varied  from  0.85  for  ^=  3  to  0.98  for  K=  30.  The  depth  de¬ 
creased  from  upstream  to  downstream  due  to  infiltration. 
The  infiltration  depth  was  zero  at  the  upstream  end  and  in¬ 
creased  to  0.5  at  the  downstream  end. 


Approximate  diffusion-wave  solution 

An  approximate  analytical  solution  of  Eq.  (25)  can  be  ob¬ 
tained  following  Parlange  et  al.  (1989).  To  that  end  Eq. 
(25)  can  be  written  as 

^  =  (1“/*-^*)^  (28) 
where  £=  1/KFq. 

Using  the  transformation  z.,  =  1  ,  Eq.  (28)  can  be 

written  as 
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□  i-ffusion  uuaue  approximation  :  Froude  number  =  0.50 

Upstream  boundary  condition  :  Constant  discharge  =  1 
Lateral  in-flom  q  =  0.00  Infiltration  F  =  0.50 

0-0  K  =  3.00  K  =  5.00  +- 1  K  =  10.00  >0<  K  =  30.00 


Fig.  2  Dimensionless  flow  depth  by  the  diffusion-wave  approxima¬ 
tion  as  a  function  of  dimensionless  distance  when  the  upstream  dis¬ 
charge  is  1.0,  Froude  number =0.5,  and  K  is  variable 


This  is  a  linear  differential  equation  and  has  an  analytical 
solution,  hti  y-h\ .  Equation  (33)  can  be  written  as 


3zi 


■z.) 


=  0. 


(34) 


f^eh^^  +  hl=zl.  (29) 

dz* 

Taking  a  Taylor-series  approximation  for  near  z*=0 

(30) 

where  HOT  represents  higher  order  terms,  and  /2»o  =  ^»  ((5)- 
From  Eq.  (29)  at  z,=0,  d/z*/(k* =-!/(/*  £).  Equation  (30) 
now  becomes 

K=K^-j^-  (31) 


Its  analytical  solution  is 


C  (/*  £  ^*0  “  ^0  +  3  z* 

-3/*e/j*oZ*  =/i^ 


where  c  is  the  constant  of  integration.  Its  value  can  be  found 
by  using  the  condition  at  z^=0 


Kozll^ 


(35) 


Substituting  Eq.  (35)  into  Eq.  (34)  yields  the  solution 


This  is  a  good  approximation  of  Eq.  (29)  near  z* = 0  but 
clearly  a  poor  approximation  when  z*  approaches  1 .  Equa¬ 
tion  (29)  can  be  written  as 


£/*/»*  d/i  3  _  _2 

— ^ —  +  «*  =  z*  .  (32) 

3  dz^ 

Substitution  of  Eq.  (31)  in  Eq.  (32)  to  replace  the  linear 
term  yields 


dhj  ^  3hi _ 3z| 

dz^  (.^  £^*0  ”  z,)  (./*  £^*q  Zj(e) 


^^^^o+3z|  (36) 

/*  ^  ^*0 

The  value  of  h^Q  is  found  by  substituting  at  z*  =  1. 

Carrying  out  the  above  substitution  and  some  algebraic  ma¬ 
nipulation  yields 

/i  £^  l4o  -  3/i  hio  +  3/*  £h^o 

+(-/2,l/*V-l)/t,o+/*2e2=0. 


(33) 


(37) 
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Dl-f-fusion  uuaue  approximation  :  Froude  number  =  0.50 

Upstream  boundary  condition  :  Constant  discharge  =  1 
Analytical  solution 

Lateral  in-flotu  q  =  0.00  Infiltration  F  =  0.50 
O©  K  =  10.00  K  =  30.00 


Fig.  3  Dimensionless  flow  depth  by  approximate  analytical  diffu¬ 
sion-wave  solution  as  a  function  of  dimensionless  distance  when  the 
upstream  discharge  is  1.0,  Froude  number =0.5,  and  AT  =  10,  30 


For  a  given  value  of  A,*  the  above  equation  must  be 
solved  to  obtain  the  value  of /t*o. 

This  is  an  explicit  solution,  for  all  quantities  on  the  right 
hand  side  are  known.  It  was  found  that  Eq.  (37)  yielded 
positive  real  roots  only  for  large  values  ofKF^.  The  solu¬ 
tion  in  terms  of  dimensionless  depth  for  two  values  of  K 
(K=  10  and  30)  is  shown  in  Fig.  3. 


Dynamic-wave  solution 


Equation  (10),  subject  to  Eq.  (15),  has  the  solution  given 
by  Eq.  (20).  Substitution  for  from  Eq.  (15)  in  Eq.  (1 1) 
leads  to 


dx„ 


--t-- 


dbc„. 


=  K 


1  — 


ui 


.  (38) 


Carrying  out  the  required  differentiation  and  simplifying 
lead  to 


KF^  ul  +  M,,(l  - X*) 


dx^ 


(i-f*xJ[Fiul  -(l-4xj] 


(39) 


Similarly,  in  terms  of  ,  Eq.  (11)  becomes 


dh,  _  K  Fj  hi  -KF^(l-f,xJ^  +  Fj 

It  is  seen  from  Eq.  (30)  that  the  depth  gradient  becomes  in¬ 
finite  when 

/'*  =  [F'o(l-/*x*)F^^  (41) 

and  is  given  by  Eq.  (27)  when  a  =/*  x*  =  1 . 

The  procedure  outlined  in  section  2  was  used  to  inte¬ 
grate  Eq.  (40)  to  get  the  value  of  for  different  values  of 
X*  for  all  12  different  cases.  For  a  sample  case  (Fo  =  0.5), 
the  solution  in  terms  of  the  dimensionless  depth  is  shown 
in  Fig.  4.  As  can  be  seen  from  the  figure,  for  a  constant 
downstream  depth  of  =  0.04,  the  upstream  depth  varied 
from  0.88  for  K=3  to  0.98  for iir=30.  The  depth  decreased 
from  upstream  to  downstream  due  to  infiltration  which  was 
zero  at  the  upstream  end  and  increased  to  0.5  at  the  down¬ 
stream  end. 


Determination  of  error 

Error  in  KW  approximation 

Equation  ( 1 8)  is  the  error  differential  equation  with  5  spec¬ 
ified  by  the  flow  depth.  For  the  KW  approximation,  is 
explicitly  given  by  Eq.  (22),  and  hp  by  the  solution  of  Eq. 
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Dynamic  uuQ^/'e  approximation  :  Froude  number  =  0,50 

Upstream  boundary  condition  :  Constant  discharge  =  1 
Lateral  intloiu  q  =  0.00  Infiltration  F  =  0.50 

K  =  3.00  K  =  5.00  +-+  K  =  10.00  >0(  K  =  30.00 


Fig,  4  Dimensionless  flow  depth  by  the  dynamic-wave  approxima¬ 
tion  as  a  function  of  dimensionless  distance  when  the  upstream  dis¬ 
charge  is  1.0,  Froude  number =0.5,  and  K  is  variable 


df  _  -2/^(£  +  l) 
ck*  3  hp 


(45) 


(40).  Substituting  for  h^,  and  in  Eq.  ( 1 8)  leads  to  the  fol¬ 

lowing  error  differential  equation: 


d£  _  (£■+!)  ,.d/iD 


(42) 


To  that  end,  dhf^/da^  is  obtained  by  differentiating  Eq. 
(22); 


2/* 


-1/3 


(43) 


and  dhp/dx^  is  given  by  Eq.  (40).  Substituting  Eqs.  (22), 
(40),  and  (43)  into  Eq.  (42)  yields  the  following  differen¬ 
tial  equation  for  error  in  the  KW  solution: 


d£  _  (£  +  1) 

(1-/.;c,)2/3 


-2/* 


3(l-/,xJ 


1/3 


(£-H)2 


2/3 


KFq^  hi  -KFq^  (1  -/,  +  £o"  ho  f,  (1  -/,  X,) 

hh-Fi{l-f,xJ^ 


.  (44) 


(£  +  l)^r££o^-/i:£o^(£  +  l)^+4£o^/ti^'^^' 
hx  [  l-Fi(E  +  lf 

where /ij^= ( 1  -/^  Equation  (45)  specifies  errorin  the 

KW  solution  as  a  function  of  distance. 

An  explicit  solution  for  error  in  the  KW  approximation 
is  not  tractable.  However,  a  numerical  solution  is  relatively 
simple.  The  longitudinal  variation  of  error  in  depth  for  the 
case  £o=0.5  is  shown  in  Fig.  5.  An  analysis  of  the  error 
profile  shows  that  the  KW  approximation  is  an  excellent 
approximation  for  large  values  of  ££o^(Ar£o‘>7.5).  The 
error  profile  is  almost  flat  for  KFq=1.5  with  error  mag¬ 
nitudes  of  1.2%.  The  error  profile  is  almost  flat  for 
KFq=1.5  with  error  magnitudes  of  1.2%.  The  error  varied 
from  1,5%  to  44%  for  KFq  =0.75.  This  indicates  that  the 
KW  approximation  is  a  poor  approximation  for  small  val¬ 
ues  of  KFq  .  It  is  also  observed  that  the  error  becomes  very 
large  at  the  downstream  end.  This  is  due  to  the  fact  that  the 
KW  approximation  require  that  the  depth  be  zero  at  the 
downstream  end,  whereas  a  small  but  finite  depth  is  re¬ 
quired  for  the  DYW  approximation. 


Error  in  DW  approximation 


Note  that  hp  is  expressed  in  terms  of  through  Eq.  (19).  The  error  differential  equation  for  DW  approximation  can 
Upon  simplifying  Eq.  (44),  be  obtained  by  substituting  Eqs.  (25)  and  (40)  in  Eq.  (42). 


in  depth 


Error  in  kinematic  approximation  :  Froude  number 

Ups  trea.m  boundary  condition  :  Constant  discharge  “  1 
Lateral  in-floiu  q  =  0,00  In-filtration  F  =  0.50 

0-0  K  =  3.00  ^  K  =  5.00  -f— h  K  =  10.00  K  =  30.00 


Dimensionless  distance  X 


Error  in  dl-Ffusion  oia'v’e  approximation 

upstream  boundary  condition  :  Constant  discharge  =  1 
Rnalytical  solution 

Lateral  in-flou  q  =  0,00  Infiltration  F  =  0.50 

0-0  K  =  10.00  K  =  30.00 


Froude  number 


limension 


less  distance  X 
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Error  In  diffusion  uudue  approximation  :  Froude  number 

upstream  boundary  condition  :  Constant  discharge  “  1 
Lateral  in-flouj  q  =  0.00  In-filtration  F  =  0.50 

0-0  K  =  3.00  A-A  K  =  5.00  -t—f-  K  =  10.00  K  =  30.00 


0.50 


C\J 


Fig.  7  Error  in  the  flow  depth  by  the  analytical  diffusion-wave  so¬ 
lution  as  a  function  of  dimensionless  distance  when  the  upstream  dis¬ 
charge  is  1.0,  Froude  number=0.5,  and  K  is  variable 


With  these  substitutions  we  obtain 


AE  _{E  +  \)\KF^hl-KF^{\-f^x^f 
dj;*  hfc  h\ 


{E  +  \f 

hK 


K F^  hl-K F^  (1  - X,  f  +  Fjhp 

hh-E^(l-f^x,f 


.(46) 


shows  that  the  DW  approximation  is  an  excellent  approx¬ 
imation  for  large  values  of  KFq  (KFq  >7.5).  The  error  pro¬ 
file  is  almost  flat  for  KFq  ^7.5  with  error  magnitudes  of 
0.8%.  The  error  varied  from  1 .42%  to  6.0%  for  KFq  =  0.75. 
As  can  be  seen  from  the  figure,  the  error  magnitude  be¬ 
comes  large  (=23%)  near  the  downstream  end  of  the  chan¬ 
nel  (1.99<.r^<2.0).  This  indicates  that  the  DW  approxi¬ 
mation  is  valid  over  a  large  range  of  KFq  and  in  most  of 
the  channel  reach. 


On  simplifying, 


d£ 

'KFihjc-KFiif^xJ' 

(£-H) 

'KFi  hi  +  Fi  hK  f,  (1  -/,  X*)  (£-H)2  -KFi(l-l  x,  f  (E+lf ' 

I 

_ 1 

hl-Fi{l-f^x,f(E  +  l)^ 

(47) 


The  longitudinal  variation  of  error  in  depth  for  the  case 
Fo=0.5  is  shown  in  Fig.  6.  An  analysis  of  the  error  profile 

Error  in  approximate  DW  solution 


Fig.  5  Error  in  the  flow  depth  by  the  kinematic-wave  approxima¬ 
tion  as  a  function  of  dimensionless  distance  when  the  upstream  dis¬ 
charge  is  1.0,  Froude  number=0.5,  and  K  is  variable 

Fig,  6  Error  in  the  flow  depth  by  the  diffusion-wave  approximation 
as  a  function  of  dimensionless  distance  when  the  upstream  discharge 
is  1.0,  Froude  number=0.5,  and  K  is  variable 


Equation  (18)  is  the  error  differential  equation  with  S  spec¬ 
ified  by  the  flow  depth.  For  the  approximate  analytical  so¬ 
lution,  h  given  by  Eq.  (36).  Substitution  for  hp  and 
in  Eq.  (18)  gives  the  differential  equation  for  error. 

To  that  end,  dhf^/dx^  is  obtained  by  differentiating  Eq. 
(36): 
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^ _ />(! 

dr*  h.  (/*  £  /i*o  +  /*  X*  -  1)  hhf^e  h^o  +  /*  x*  - 1) 

and  d/2£)/dr*  is  given  by  Eq.  (39). 

Substituting  Eqs.  (36),  (40),  and  (48)  into  Eq.  (42)  yields  the  following  differential  equation  for  error  in  the  approx¬ 
imate  analytical  solution. 


d£  _  (£  + 1) 

(£  +  1)' 

K Fi hl-KF^{l-f.x.f{E  + 1)4  foV.  hK (1  -f.x.){E  +  lf 

dj*.  hi 

f.£h.o  +  f.x.-l) 

hK 

[hl-F^(l-f.x.)HE  +  lf] 

where  hf^  is  given  by  Eq.  (36).  The  longitudinal  variation 
of  error  is  shown  in  Fig.  7.  The  magnitude  of  error  varied 
from  about  6.4%  at  the  upstream  end  to  greater  than  30% 
at  the  downstream  end.  This  shows  that  the  approximate 
analytical  diffusion-wave  solution  is  not  suitable. 


Conclusions 

The  following  conclusions  can  be  drawn  from  this  study: 
( 1 )  When  the  upstream  condition  of  finite  discharge  and 
depth  is  known,  the  governing  differential  equation  can¬ 
not  be  integrated  directly  from  upstream  to  downstream. 
An  iterative  procedure  must  be  followed  by  which  the 
downstream  depth,  K  and  Fq  are  varied  until  the  upstream 
depth  agrees  with  the  given  condition.  (2)  For  upstream 
boundary  conditions  of  finite  discharge  and  depth,  the  KW 
approximation  was  in  excellent  agreement  with  the  DYW 
representation  for  large  values  of  KFq.  The  error  profiles 
were  flat  [almost  a  constant  value  of  error  (=1.2%)]  for 
KFq  =7.5  and  curved  downwards  for  smaller  values  of  KFq 
with  errors  varying  between  18%  and  44%  for  ^Fo^=0.75. 
The  KW  approximation  is  a  good  approximation  for  large 
values  of  KFq,  but  a  poor  approximation  for  small  values 
of  KFq.  (3)  For  upstream  boundary  conditions  of  finite 
discharge  and  depth,  the  DW  approximation  was  in  excel¬ 
lent  agreement  with  the  DYW  representation  for  large  val¬ 
ues  of  KFq  in  the  region  0.0<jc^<  1.98.  The  error  profiles 
were  flat  [almost  a  constant  value  of  error  (=0.8%)]  for 
KF^  =  1.5  and  =1.4%  for  KF^  =  0.75.  The  DW  approxima¬ 
tion  was  found  to  be  a  good  approximation  for  the  range 
of  values  of  KFq  considered  (0.7 5  < KFq  <7. 5)  and  in  the 
region  0.0<x*<1.98.  (4)  The  approximate  analytical  so¬ 
lution  performed  poorly  for  larger  values  of  KFq  and  a  vi¬ 


able  solution  could  not  be  obtained  for  small  values  of 

Abbreviations:  C=^hezy’s  resistance  parameter  f^rate  of 

infiltration  (L/T);  f  =  average  rate  of  infiltration  {L/T)\  F=Froude 
number;  Fo= normalizing  Froude  number;  F*  =  Froude  number  in  di¬ 
mensionless  domain;  acceleration  due  to  gravity  (L/T^);  /z  =  fIow 
depth  (L);  /i*= dimensionless  flow  depth;  //o= norma!  depth  of  flow 
(L);  A’=  kinematic  wave  number;  L= length  of  the  channel  (L);  Q- 
discharge  per  unit  width  (L^I(TL))\  Qo  =  *^ormal  discharge  (I}i{tL))\ 
2*  =  dimensional  discharge;  frictional  slope;  5o=bed  slope; 

r=  duration  of  irrigation  (T)\  normalizing  time  (T)\  M  =  flow  ve¬ 

locity  (L/T);  t/o= normal  velocity  {UT)\  dimensionless  veloc¬ 
ity;  .v= distance  measured  from  the  upstream  boundary  (L);  jCo= nor¬ 
malizing  distance  (L);jc*= dimensionless  distance;  /5=  resistance  pa¬ 
rameter  (T^/L). 
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ABSTRACT:  Hydrodynamic  models  employed  for  computing  flood  discharges  are  based  on 
the  shallow  water-wave  theory  that  is  described  by  the  St.  Venant  (SV)  equations. 
These  models  are  derived  from  either  the  kinematic— wave  (KW)  approximation,  the 
dif fusion- wave  (DW)  approximation  or  the  dynamic-wave  (DYW)  representation  of  the 
SW  equations.  In  the  studies  reported  to  date,  different  types  of  criteria  have 
been  established  to  evaluate  the  adequacy  of  the  KW  and  DW  approximations,  but  no 
explicit  relations  either  in  time  or  in  space  between  these  criteria  and  the  errors 
resulting  from  these  approximations  have  been  derived  yet.  Furthermore,  when  doing 
hydrologic  modeling,  it  is  not  evident  if  the  KW  and  DW  approximations  are  valid  on 
one  hand  for  the  entire  hydrograph  or  a  portion  thereof,  and  on  the  other  hand  for 
the  entire  channel  length  or  a  portion  thereof.  In  other  words,  most  of  these 
criteria  take  on  fixed  point-values  for  a  given  rainfall-runoff  event.  This  paper 
attempts  to  derive,  under  simplified  conditions,  error  equations  for  the  KW  and  DW 
approximations  for  space- independent  as  well  as  for  time— independent  flows,  which 
provide  a  continuous  description  of  error  in  the  flow-discharge  hydrograph.  For 
space- independent  flows,  a  dimensionless  parameter  y  is  defined  which  reflects  the 
effect  of  initial  depth  of  flow,  channel-bed  slope,  lateral  inflow,  and  channel 
roughness.  For  time-independent  flows,  the  dimensionless  parameter  is  the  product 
of  the  kinematic  wave  number  and  the  square  of  the  Froude  number.  The  kinematic 
wave,  diffusion  wave  and  dynamic  wave  solutions  are  parameterized  through  these 
parameters.  By  comparing  the  kinematic  wave  and  diffusion  wave  solutions  with  the 
dynamic  wave  solution,  equations  are  derived  in  terms  of  these  parameters  for  the 
error  in  the  kinematic  wave  and  diffusion  wave  approximations. 


1  INTRODUCTION 

Physically-based  models  of  overland 
flow,  channel  flow,  surface  irrigation, 
and  many  other  phenomena  involving 
unsteady,  free  surface  open  channel 
flows  are  based  on  the  shallow  water 
wave  (SWW)  theory.  These  models  are 
based  either  on  the  kinematic  wave  (KW) 
approximation  (Lighthill  and  Whitham, 
1955) ,  diffusion  wave  approximation 
(DW) ,  or  dynamic  wave  (DYW)  representa¬ 
tion.  Lighthill  and  Whitham  (1955) 
showed  that  at  the  Froude  numbers  less 
than  one  (appropriate  to  flood  waves) 
the  dynamic  waves  are  rapidly  atten¬ 
uated  and  the  kinematic  waves  become 
dominant.  Using  a  dimensionless  form 
of  the  St.  Venant  (SV)  equations, 
Woolhiser  and  Liggett  (1967)  obtained 
what  is  now  referred  to  as  the  kinema¬ 
tic  wave  number,  K,  as  a  criterion  for 
evaluating  the  adequacy  of  the  KW 
approximation.  For  K  greater  than  20, 
the  KW  approximation  was  considered  to 
be  an  accurate  representation  of  the  SV 


equations  in  modeling  of  overland  flow. 
However,  no  relation  between  K  and  the 
error  in  the  KW  approximation  was  sug¬ 
gested.  Morris  and  Woolhiser  (1980) 
modified  the  above  criterion  with  an 
explicit  inclusion  of  Froude  number, 

Fq,  and  showed,  based  on  numerical 
experimentation,  that  Fq^K  >  5  was  a 
better  indicator  of  the  adequacy  of  the 
KW  approximation.  A  relation  between 
this  criterion  and  the  error  resulting 
from  the  KW  approximation  was  not 
derived,  however. 

Using  a  linear  perturbation  analy¬ 
sis,  Ponce  and  Simons  (1977)  derived 
properties  of  the  KW  and  DW  approxima¬ 
tions  as  well  as  DYW  representations  in 
modeling  of  open  channel  flows.  They 
derived  a  spectrum  showing  the  regions 
of  the  validity  of  the  KW  and  DW 
approximations.  Menendez  and  Norscini 
(1982)  extended  the  work  of  Ponce  and 
Simons  by  including  the  phase  lag 
between  the  depth  and  velocity  of  flow. 
Their  results  were,  however,  similar  to 
those  of  Ponce  and  Simons  (1977) .  In 
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another  but  similar  study,  Ponce,  et 
ai.  (1978),  based  on  propagation  char¬ 
acteristics  of  sinusoidal  perturbation, 
derived  criteria  to  evaluate  the  ade¬ 
quacy  of  the  KW  and  DW  approximations. 
Daluz  Viera  (1983)  compared  solutions 
of  the  SV  equations  with  those  of  the 
KW  and  DW  approximations  for  a  range  of 
Fq  and  K,  and  defined  the  regions  of 
validity  of  these  approximations  in  the 
K-Fq  space. 

Fread  (1985)  developed  criteria  for 
defining  the  range  of  application  of 
the  KW  and  DW  approximations.  These 
were  based  on  an  analysis  of  the  magni¬ 
tude  of  the  normalized  errors  in  the 
moment-urn  equation  due  to  omission  of 
certain  terms.  In  a  comprehensive 
study,  Ferrick  (1985)  defined  a  group 
of  dimensionless  scaling  parameters  to 
establish  the  spectrum  of  river  waves, 
with  continuous  transitions  between 
wave  types  and  subtypes.  With  the  aid 
of  these  parameters  he  was  able  to 
discern  when  the  KW  and  DW  approxima¬ 
tions  would  be  valid. 

In  most  of  these  studies,  different 
types  of  criteria  have  clearly  been 
established  to  evaluate  the  adequacy  of 
the  KW  and/or  DW  approximations,  but  no 
explicit  relations  either  in  time  or 
space  between  these  criteria  and  the 
errors  resulting  from  these  approxima¬ 
tions  have  been  derived  yet.  Further¬ 
more,  when  doing  hydrologic  modeling  it 
is  not  evident  if  the  KW  and  DW  approx¬ 
imations  are  valid  for  the  entire  flood 
hydrograph  or  a  portion  thereof.  In 
other  words,  most  of  these  criteria 
take  on  fixed  point  values  for  a  given 
event.  The  objective  of  this  study  is 
to  derive,  under  simplified  conditions, 
error  equations  for  the  KW  and  DW 
approximations  for  space- independent  as 
well  as  time- independent  flows,  which 
specify  errors  as  a  function  of  time. 


2  SHALLOW  WATER-WAVE  (SWW)  THEORY 


rainfall  intensity  (L/T) ,  f  is  uniform 
infiltration  rate  (L/T) ,  g  is  accelera¬ 
tion  due  to  gravity,  x  is  space  coordi¬ 
nate  in  the  direction  of  flow  (L) ,  t  is 
time  (T) ,  Sq  is  bed  slope,  and  is 
frictional  slope.  Note  Q  «  uh  is  dis¬ 
charge  (L^/TL)  per  unit  width.  can 

be  approximated  as 


(3) 


where  p  is  some  resistance  parameter. 

If  the  Chezy  relation  is  used  for 
representing  the  friction  then  P  » 
g/C^,  where  C  is  Chezy' s  resistance 
parameter. 

The  DYW  representation  employs  the 
full  form  of  equations  (1)  and  (2) , 

The  KW  approximation  is  based  on  equa¬ 
tion  (1)  and  equation  (2)  with  the  left 
side  omitted, 

-  V  -  ?  =  ° 


The  DW  approximation  uses  equation  (1) 
and  equation  (2)  with  local  and  convec¬ 
tive  acceleration  deleted. 


9(S„  -  S,) 


qu 

h 


(5) 


Analytical  solutions  of  the  SV 
equations  or  their  variants  in  the  KW 
and  DW  approximations  are  tractable 
only  for  simple  cases.  To  that  end, 
both  space-independent  and  time- 
independent  cases  are  considered  in 
this  study.  In  the  first  case,  the 
water  surface  is  flat. 


3  SPACE- INDEPENDENT  FLOWS 

3.1  Governing  Equations 

For  space- independent  (or  uniform) 
flows,  equation  (1)  takes  the  form 


The  SWW  theory  can  be  described  by  some 
form  of  the  SV  equations.  For  flow 
over  an  infiltrating  plane  subject  to 
uniform  rainfall,  these  equations  can 
be  written  in  one-dimensional  form  on  a 
unit  width'  basis  as 


continuity  equation, 
q  -  f 


(1) 


momentum  equation. 


s  'I  ■  ‘^'V  =f>  -  ¥ 

(2) 

where  h  is  the  depth  of  flow  (L) ,  u  is 
local  mean  velocity  (L/T) ,  q  is  uniform 


^  = 
dt 


q  “ 


f 


(6) 


and  equation  (2)  becomes 


^  =  g(S  -  S  ) 
dt  0  f 


qu 

h 


(7) 


Equations  (6)  and  (7)  are  the  governing 
equations  for  the  DYW  representation 
for  spatially  uniform  flows.  The  KW 
approximation  is  based  on  equation  (6) 
and  equation  (7)  with  the  left  side 
dropped. 


The  DW  approximation  uses  equation  (6) 
as  well  as  equation  (8) .  Therefore, 
for  spatially  uniform  flows,  the  KW 
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(14) 


approximation  is  identical  to  the  DW 
approximation.  Equation  (8)  can  also 
be  approximated  by  neglecting  the 
momentum  exchange  between  lateral 
inflow  and  longitudinal  channel  flow  as 

So  *  Sf  (9) 


and  equation  (2)  becomes 

<iu2  +  gh)  =  g(SQ  -  S^)  -^(15) 


which  can  be  expressed  as  equation  (3) . 
Similarly,  equation  (7)  can  be  written 
as 


du  =  g(S  - 
dt  ^  0 


(10) 


Depending  upon  the  presence  or 
absence  of  f,  equation  (6)  can  also  be 
simplified.  If  f  =  0,  then 

#  -  q  (11) 


3.2  Types  of  Scenarios 

Depending  upon  the  presence  of  lateral 
inflow  and  infiltration,  four  different 
scenarios  can  be  considered:  (1)  f  =» 

Of  ^  ^  constant;  this  includes  the 

case  q-0.  (2)  q=qo=  constant,  and 

f  ~  f0  =  constant;  and  this  includes 
the  case  q  -  f  -  0.  (3)  q  -  f  =  0,  q  - 

qo  =  constant;  this  includes  the  case 
q  =  f  -  0.  (4)  q  =  0,  f  =  fo  =  con~ 
stant;  this  includes  the  case  f  =  0. 

It  may  be  noted  that  the  scenario  with 
q  =  0  in  equation  (11)  or  (q  -  f)  =0 
in  equation  (6)  applies  to  the 
recession  hydrograph.  The  same  applies 
if  (q  -  f)  <0. 


Equations  (14)  and  (15)  are  the  govern¬ 
ing  equations  for  the  DYW  representa¬ 
tion  for  time- independent  flows.  The 
KW  approximation  is  based  on  equation 
(14)  and  equation  (15) ,  with  the  left 
side  deleted,  leading  to  equation  (8) . 
The  DW  approximation  uses  equation  (14) 
and  equation  (15)  with  convective 
acceleration  deleted. 


qu 

h 


(16) 


Equations  (15)  and  (16)  can  also  be 
approximated  by  neglecting  the  momen¬ 
tum-exchange  between  lateral  inflow  and 
longitudinal  channel  flow  respectively 
as 


^  (lu2  .  ,h)  .  5(S„  -  S^l  <17, 


and 


dh  =  c; 
dx  0 


(18) 


Depending  upon  the  presence  of  f,  ecjua- 
tion  (14)  can  also  be  simplified.  If 
f  -  0 ,  then 

==  q  (19) 


3.3  Initial  Conditions 


Two  types  of  initial  conditions  (t  =  0) 
can  be  assumed:  , 


(1)  h(0)  =  hQr  u(0)  -  uq,  (12) 

(2)  u(0)  =  0,  h(0)  -  0  (13) 


3.4  Scenarios  for  Determination  of 
Error 


4.2  Types  of  Scenarios 

Depending  upon  the  presence  of  lateral 
inflow  and  infiltration,  four  different 
scenarios  can  be  considered  as  in  the 
preceding  section.  It  may  be  noted 
that  the  scenario  with  q  =  0  in  equa¬ 
tion  (19)  or  q  -  f  =  0  in  equation  (14) 
applies  to  the  case  of  losing  flow. 

This  same  would  apply  if  q  -  f  <  0 . 


Error  equations  have  been  derived  by 
Singh  (1992a,  1992b)  for  the  KW  and  DW 
approximations  under  the  above-men¬ 
tioned  conditions  for  four  different 
scenarios .  Treated  together,  eighteen 
cases  result  and  are  summarized  in 
Table  1. 


4  TIME- INDEPENDENT  FLOWS 

4.1  Governing  Equations 

For  time-independent  flow,  equation  (1) 
takes  the  form 


4 . 3  Boundary  Conditions 

Two  types  of  conditions  at  the  upstream 
boundary  (x  =  0)  can  be  assumed: 

(1)  h(0)  =  hQ,  u(0)  -  Uq,  (20) 

(2)  u(0)  =  0,  h(0)  =  ho  (21) 

Depending  upon  the  type  of  flow,  the 
boundary  conditions  may  have  to  be 
specified  at  the  upstream  boundary  as 
well  as  the  downstream  boundary  or  at 
the  upstream  boundary  alone. 

The  upstream  boundary  condition, 
given  by  equation  (21) ,  influences  both 
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Table  1.  List  of  cases  for  error  analysis  for  space-independent  flows 
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the  supercritical  and  stibcritical  flow 
outside  zone  A  of  the  characteristic 
solution  domain.  The  downstream  boun¬ 
dary  condition  (x  »  L)  can  be  of  two 
types: 


6.1  Kinematic  Wave  and  Diffusion  Wave 
Solution 

Equation  (11),  subject  to  equation 
(12),  has  the  solution: 


(1)  Critical  flow  downstream 
boundary  condition: 


=  hg  +  qQt 


(27) 


u(L)  -  [gh(L)]0-5 

where  L  is  the  channel  length.  This 
occurs  when  the  channel  ends  at  the 
steep  bank  of  a  river.  For  supercri¬ 
tical  slow 


u(L)  >  [h(L)g]0-5  (23) 

(2)  Zero-depth  gradient  downstream 
boundary  condition: 


dhllil  =  0 
dx 


(24) 


4.4  Scenarios  for  Determination  of 
Error 


From  the  kinematic  wave  approximation. 


u 


0.5 


h 


0.5 


(28) 


It  is  convenient  to  define  a  dimension¬ 
less  time  T  as 


ho  +  qQt 
^0  ' 


0, 


X  >  1  (29) 


Elation  (14)  can  be  expressed  in 
dimensionless  form  as 


E<^ation  (34)  can  be  expressed  in 
dimensionless  form  as 


Error  equations  have  been  derived  by 
Singh  and  Aravamuthan  (1992a,  1992b, 
1992c)  for  the  KW  and  DW  approximations 
under  the  above-mentioned  conditions 
for  four  different  scenarios.  Taken 
together,  these  will  give  rise  to  21 
cases  that  are  siammarized  in  Table  2. 


v-§=t0-5,  u-(£|!10)0-5  (30) 

The  kinematic  wave  (KW)  or  diffusion 
wave  (DW)  solution  is  given  by  equa¬ 
tions  (33)  and  (34).  The  velocity 
increases  parabolically  with  T. 


5  DEFINITION  OF  ERROR 


The  relative  error  E  is  defined  as 


(25) 


where  Sj^  is  the  solution  from  the  KW  or 
DW  approximation,  and  Sq  is  the  solu¬ 
tion  from  the  DYW  representation.  The 
solution  can  be  either  in  terms  of 
depth  (h),  velocity  (u) ,  or  discharge 
(Q) ,  i.e.,  S  {u,h,Q).  Subscripts  K  and 
D  correspond  to  the  KW  (or  DW)  and  DYW 
solutions,  respectively.  The  differen¬ 
tial  equation  of  E  can  be  obtained  by 
differentiating  equation  (25)  as 

^  ,  iE+il  ^  _  jEil)2  dSo 

dx  Sk  dx  Sk  dir 

The  differential  equation  for  error 
can,  however,  be  defrned  without  expli¬ 
citly  knowing  Sn,  so  long  as  Sv  is 
explicitly  known.  ^ 


6  ERROR  EQUATIONS  FOR  SPACE- 
INDEPENDENT  FLOWS 

We  consider  the  case  where  f  -  0,  and 
^  ^  ^0- 


6.2  Dynamic  Wave  Solution 

Equation  (11)  has  the  solution  given  by 
equation  (27) .  Equation  (2)  can  be 
expressed  in  dimensionless  terms  as 

dv  ^  _  yO.5  ^  ^  ^  4g^(JSQhQ 

dt  2  2  X  '  ^  *  2 

^0 


(31) 

Equation  (31)  is  a  special  case  of 
the  Riccati  equation.  When  t  =*  1, 

V  =»  1,  and  dv/dx  -  0.  With  these 
initial  conditions  in  hand,  equation 
(31)  was  solved  by  using  the  4th  order 
Runge— Kutta  method.  For  a  fixed  y,  v 
increases  with  increasing  x,  and  for  a 
fixed  X,  it  increases  with  y.  However, 
for  y  ^  1.5,  V  is  not  very  sensitive  to 


6.  3  Error  in  KW  and  DW  Approximations 

By  making  use  of  equations  (30)  and 
(31)  in  eqfuation  (26),  one  obtains  the 
error  differential  equation: 


+  C2(Y,T)  e2 


(32a) 
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Table  2.  List  of  cases  for  error  analysis  for  time- independent  flows. 
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(17)  I  Eqs.  (14)  and  (17)  |  Eg.  (20)  |  N.A. 


Diffusion  wave  approximation 


gtrO’O  9r0'0-  Z.  1 1*0-  S6 1  *0-  6Z.rO-  09r0-  235*0-  509*0- 189*0-  594*0- 


pD9q  ui  jojja  96d:^u0OJ9^ 
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Dimensionless  distance  X 

Figure  2.  Error  in  the  KW  approximation  as  a  function  of  dimensionless  distance  for  time 
independent  flows. 


Kinematic  wave  approximation  :  Froude  number 

Upstream  boundary  conditions:  Discharge  M  =  0.00  Depth  H  =  1 .03 

Lateral  inflow  q  =  1 .00  Infiltration  F  =  0.00 


Dimensionless  distance 


then 


h 


* 


h  ,  qpx  +  a  2/3 

H  Uq^’O 


(39) 


7.2  Diffusion  Wave  Solution 


Equation  (1)  has  the  solution  given  by 
equation  (37) .  Equation  (2)  reduces  to 


(40) 


9  ERROR  EQUATIONS  FOR  TIME- INDEPENDENT 
FLOWS:  CRITICAL  FLOW-DEPTH  DOWN¬ 

STREAM  BOUNDARY  CONDITION 

The  error  equations  for  the  critical 
depth  downstream  boundary  condition 
were  obtained  following  the  procedure 
of  Section  7.  For  the  DW  and  KW 
approximations r  the  errors  obtained 
were  almost  of  the  same  magnitude  as  of 
those  obtained  using  the  zero  depth 
gradient  downstream  condition  in  the 
region  0.1  <  x  <  1.0. 


With  the  upstream  boundary  condition 
given  by  equation  (20) ,  equation  (40) 
can  be  written  in  terms  of  h  as 


^  =  s  - 

dx  0 


(41) 


7.3  Dynamic  Wave  Solution 

Equation  (1)  has  the  solution  given  by 
equation  (44) .  Equation  (2)  reduces  to 

^  (1  »  ghi  -  -  gP  ^  (42) 

which  can  be  expressed  in  terms  of  h  as 


dh  ^  _ 1 _ 

^  gh^  -  (q^x+a)^ 


[gSQh3 


-  qp  “  P(qpx+a)2] 


(43) 


7.4  Errors  of  KW  and  DW  Approximations 

The  DW  approximation  gave  good  results 
for  the  entire  channel  reach  with  error 
less  than  1%.  The  err, or  decreased  with 
increasing  values  of  KFq^,..  as  shown  in 
Figure  2.  The  KW  approximation  gave 
adequate  results  (error  <  10%I  in  the 
region  (0.1  <  x  <  1.0)  for  Fq^  values 
greater  than  30,  as  shown  in  Figure  3. 


10  CONCLUSIONS 

For  space-independent  flows,  the  kine¬ 
matic-wave  and  dif fusion- wave  approxi¬ 
mations  are  sufficiently  accurate  when 
the  dimensionless  parameter  y  ^  3. 

This  parameter  reflects  the  effect  of 
initial  flow  depth,  bed  slope,  lateral 
inflow,  and  channel  roughness.  The 
error  of  these  approximations  declines 
exponentially  when  the  dimensionless 
time  exceeds  5.  The  dimensionless  time 
is  obtained  with  the  use  of  initial 
depth . 

For  time- independent  flows  with  a 
constant  depth  upstream  boundary  con¬ 
dition  and  with  parameters  (f  -  0.0, 
a  ^  0.0,  H  »  1.03,  Fq  -  1.0),  the  DW 
approximation  was  quite  accurate,  with 
errors  in  depths  of  less  than  1%,  The 
errors  decreased  with  an  increase  in 
the  value  of  KFq^.  The  KW  approxima¬ 
tion  was  accurate  with  errors  of  less 
than  10%  in  the  region  (0.1  <  x  <  1.0) 
and  for  KFq^  values  greater  than  30. 

For  smaller  values  of  KFq^,  the  KW 
approximation  was  not  adequate. 
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8  ERROR  EQUATIONS  FOR  TIME- INDEPENDENT 
FLOWS:  ZERO  DEPTH-GRADIENT  DOWN¬ 

STREAM  iBOUNDARY 

Following  the  procedure  of  Section  7, 
error  equations  were  derived  for  the 
case  of  the  zero-depth  gradient  down¬ 
stream  boundary  condition.  The  DW 
approximations  gave  adequate  results  in 
the  region  (0.1  <  x  <  1.0)  and  for  KFq^ 
values  greater  than  30.  For  larger 
values  of  KFq^,  the  errors  were  less 
than  1%  in  the  region  (0.1  <  x  <  1,0) 
and  increased  to  about  21%  at  the 
upstream  end  of  the  channel.  The  KW 
approximation  gave  reasonable  results 
(errors  <  10.1)  for  KFq^  greater  than 
30  and  in  the  region  0.1  <  x  <  1,0. 
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ABSTRACT 

Wang,  G.-T.,  and  Singh,  V.P.,  1992.  Muskingum  method  with  variable  parameters  for  flood  routing  in 
channels.  J.  Hydrol.,  134:  57-76. 

Three  versions  of  the  Muskingum  method  with  variable  parameters  for  flood  routing  in  open  channels 
were  derived.  The  reach  travel  time  was  obtained  from  simplification  of  the  St.  Venant  equations.  The  three 
versions  were  applied  to  three  example  data  sets  and  the  results  obtained  were  found  to  be  more  accurate 
than  those  of  the  conventional  Muskingum  method  with  parameters  obtained  by  trial  and  error  or 
optimization. 


INTRODUCTION 

Since  its  development  around  1934  by  McCarthy  (1938),  the  Muskingum 
method  for  routing  flood  waves  in  rivers  and  channels  has  been  a  subject  of 
many  investigations  (Singh,  1988).  The  method  represents  a  river  reach  as  a 
linear  time-invariant  system  with  its  inflow  /,  outflow  Q,  and  storage  w. 
related  as 

n’  =  K[XI+{\-X)Q]  (1) 

in  which  K  and  X  are  parameters. 

The  Muskingum  method,  based  on  eqn.  (1)  and  the  water  balance  equation, 
can  be  expressed  as 

Q{i+\)  =  C,/(/+l)-fC/(O  +  C30(/),  /•  =  1,2,...  (2) 

where 


'  Visiting  Scholar.  Department  of  Civil  Engineering.  Louisiana  State  University.  Baton  Rouge. 
LA  70803.  USA. 
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p 

> 

1 

(3) 

0.5A/+(1-A:)is: 

o.5At+ Jir/: 

0.5At  +  (l-AO/: 

(4) 

-0.5Ar-h(l-J!0^ 

0.5A/  +  (1-AO/: 

(5) 

where  At  is  the  time  interval  (/,+  ,  -/,)  or  [(/+ !)-(/)],  Qii+  1)  is  the  outflow 
at  the  end  of  the  time  interval,  /(/+ 1)  is  the  inflow  at  the  end  of  the  time 
interval,  Q{i)  is  the  outflow  at  the  beginning  of  the  time  interval,  and  /(/)  is 
the  inflow  at  the  beginning  of  the  time  interval. 

In  the  conventional  Muskingum  method,  the  parameters  K  and  X  are 
considered  constant  and  are  determined  by  calibration  using  measured  inflow 
and  outflow  hydrographs.  This  assumption  makes  the  parameters  dependent 
on  the  inflow-outflow  values  used  to  evaluate  them  (Dooge,  1973).  As  a  result, 
the  parameter  values  change  from  one  set  of  inflow-outflow  hydrograph  data 
to  another.  A  more  physically  realistic  approach  is  to  consider  the  parameters 
K  and  X  to  vary  in  time  and  space  according  to  the  flow  variability  (Miller  and 
Cunge,  1975).  Weinmann  and  Laurenson  (1979)  reviewed  approximate  flood 
routing  methods,  with  focus  on  the  Muskingum  diffusion  method.  They 
proposed  an  inequality  for  use  in  routing  calculations  if  negative  outflows 
were  to  be  ignored. 

Cunge  (1969)  derived  eqn.  (2)  using  a  finite  difference  approximation  of  the 
St.  Venant  equations  with  the  inertia  term  neglected,  and  expressed  K  and  X 
in  terms  of  physical  river  characteristics 


A.r  AxA 


(6) 

(7) 


where  A.v  is  the  reach  length  (space  interval),  c  is  the  flood  wave  celerity,  q  = 
QIB  is  the  unit  width  discharge,  Q  is  the  discharge,  B  is  the  top  width.  5o  is 
the  channel  bed  slope,  and  P  is  an  exponent. 

Koussis  (1976.  1978,  1980,  1983)  considered  the  macro  form  of  the 
continuity  equation  with  no  lateral  inflow  or  lateral  outflow.  By  carrying  out 
the  discretization  in  space  only,  retaining  continuous  functions  for  the  time 
derivatives,  using  a  weighting  scheme  for  the  time  derivative,  and  assuming 
linear  variation  of  flow  over  the  time  interval  At.  he  derived  a  more  general 
form  of  the  Muskingum  method  wherein  the  parameters  were  discharge- 
dependent.  Koussis  was  thus  able  to  account  for  the  nonlinear  nature  of  the 
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flood  wave  propagation  through  an  explicit  consideration  of  the  unsteady- 
flow  loop-shaped  storage-discharge  relationship.  Koussis  (1976)  was  among 
the  first  to  propose  a  criterion  for  space  interval  or  grid  size  to  be  used  in  the 
Muskingum  flood  routing.  However,  he  assumed  X  to  be  constant  on  the 
grounds  that  computations  were  relatively  insensitive  to  this  parameter. 

In  the  development  of  its  HEC-1  computer  program,  the  Hydrologic 
Engineering  Center  (1982)  took  A'  =  0.3  and  estimated  K  from  the  mean 
velocity.  In  order  to  make  the  routed  and  observed  hydrographs  comparable, 
K  was  adjusted  as  necessary.  Ponce  (1979)  considered  a  simplified  form  of  the 
Muskingum  method  by  taking  A'  =  0  and  K/Ai  =  1,  where  At  is  the  time 
interval.  Ponce  and  Yevjevich  (1978)  allowed  the  parameters  K  and  X  to  vary 
in  space  and  time  and  computed  them  for  each  computational  cell  using 
three-point  and  four-point  numerical  methods.  Ponce  and  Theurer  (1982) 
calculated  the  Muskingum  parameters  based  on  channel  and  grid  character¬ 
istics.  Based  on  numerical  experimentation,  they  found  an  upper  limit  to  the 
space  step  if  accuracy  was  to  be  preserved. 

In  this  study,  the  reach  travel  time  K  and  the  weighting  factor  X  were 
derived  from  a  simplified  form  of  the  St.  Venant  equations.  K  was  a  function 
of  storage,  discharge  and  other  hydraulic  characteristics,  that  is,  K  =  w/ 
^  ^  function  of  discharge  and  other  hydraulic 

characteristics.  The  methods  of  parameter  estimation  employed  here  are 
different  from  the  existing  methods.  Three  versions  of  the  Muskingum  method 
with  variable  parameters  are  suggested,  in  which  the  parameters  were 
determined  either  by  using  a  direct  search  method  for  inflow-outflow  data,  or 
by  using  physical  characteristics  of  the  river  reach.  These  versions  were  tested 
using  three  data  sets  and  were  also  compared  with  the  conventional 
Muskingum  method. 

THEORETICAL  BASIS  FOR  VARIABLE  PARAMETERS 


Unsteady  flow  in  open  channels  can  be  represented  by  the  St.  Venant 
equations  of  continuity 


^+^  =  0 
Cl  cx 


(8) 


and  momentum 

cZ  1  cV  VdV  V- 

- r—  — - r-  H - - ^  Z'’  6 

CX  g  ct  g  CX  C  R 

where  A  is  the  cross-sectional  area  of  flow.  V  is  the  average  velocity  of  the  flow 
cross-section.  Q  is  the  discharge.  Q  =  AV,  Z  is  the  water  level,  g  is  the 
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gravitational  acceleration,  R  is  the  hydraulic  radius,  and  C  is  the  Chezy’s 
coefficient.  For  many  flow  situations,  the  inertia  and  acceleration  terms  can 
be  considered  much  smaller  than  friction  slope  (Henderson,  1967),  and  can, 
therefore,  be  neglected.  Equation  (9)  then  reduces  to 


8Z  _ 

dx  ~  C^R 

(10) 

or 

V  =  C^(RS) 

(11) 

where  S  is  the  water  surface  slope. 
Average  reach  travel  time 


For  shallow  wide  rectangular  channels,  R^h,  and  eqn.  (11)  becomes 
V  =  =  c^/m  (12) 

A  =  PR^Ph 

where  P  is  the  wetted  perimeter  of  the  flow  cross  section,  and  can  be  approxi¬ 
mated  for  wide  concave  channels  as 

P  =  ah""'  (13) 

where  a  is  constant,  and  m  is  exponent  with  values  between  1.2  and  2.0. 
Therefore 

A  =  ah""'  x/i  =  ah'"  (14) 

Multiplying  both  sides  of  eqn.  (14)  by  the  reach  length.  I.  and  taking  the 
average  values  of  A  and  h  over  L  2iS  A  and  /?, 

vr  =  LA  =  La{hy"  (15) 

By  substituting  eqn.  (12)  into  eqn.  (15)  and  eliminating  /7,  the  following  is 
obtained; 


KaV-'"* 

(C-5)'” 


(16) 


According  to  hydraulics  of  open  channel  flow  (Henderson.  1967) 

Q  =  PR'-’QS„  +  SJ=  ,|7a, 

K  =  «-(S„  +  S,)'=  (17b) 

where  5„  is  the  bottom  slope  of  the  river,  =  S-5„  and  a  is  an  exponent, 
X  =  0,5  for  the  Chezy  formula,  x  -  2/3  for  the  Manning  formula. 
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By  substituting  eqn.  (17b)  into  eqn.  (17a)  and  eliminating  R,  the  following 
is  obtained: 

By  substituting  eqn.  (18)  into  eqn.  (16)  and  simplifying, 

Kn  Kn  V  n2m+\ 


+a)  /^1/(1  +a) 


C''<'+“>(5o  +  5a)' 


[(1/2(1+ a)]  /na/(l  +2) 


Ka 

^  px{lm^\)j{\  +:l) 


-2m3t);'(l +a)^^^  ^  ^  y2m 


+  1  )/[2(  I  +  a)]  /^a(2m  +  1 )/( 1  +  a) 


(19a) 


(2w+  l);[2(I  +  x)] 


2m+  l)/(l  +a) 


Then  eqn.  (1 9a)  reduces  to 

H-  = 


(1 9b) 


)x(2m+  1 )  { I  ^x) 


When  [a(2m +!)]/(!  +  a)  =  1.0,  eqn.  (20)  reduces  to  eqn.  (6)  presented  by 
Ponce  (1979). 


Dimensionless  weighting  coefficient  X 


For  the  Muskingum  method,  X  can  be  represented  (Kalinin  and  Milyukov, 
1957)  as 


where  /  is  the  characteristic  length  which  can  be  expressed  as 

/  =  =  gop(^-7/n)~ 

■^0  (^  Qi) )  *^0  -  ^Q{) 


where  //„  is  the  bottom  elevation.  H  is  the  water  level,  and  Q^^  is  the  steady- 
flow  discharge. 

For  steady  flow,  the  relationship  between  the  water  level  and  the  discharge 
is 


//-//„=/;  =  b^,Qi  (23) 

where  h  is  the  water  depth,  is  the  constant,  and  [i  is  an  exponent.  Taking 
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the  derivative  of  eqn.  (23) 

=  wt'  (24) 

oQo  cQo 

and  substituting  into  eqn.  (22)  yield 

Qq  d{H —H)  _  Qo.  cffi  Pf'i 

'  —  t: - -  —  ~^f^oPQo  —  ~^Qo  (^j) 

Oq  O^q  Oq  Oq 

Substitution  of  eqn.  (25)  into  eqn.  (21)  yields 
y  *  ^  n  <:  ^oPQo 

^=2-2l  =“-^-lZsr 

The  relationship  between  steady-flow  discharge  Qq  and  unsteady-flow 
discharge  Q  can  be  expressed  (Henderson,  1967)  as 

o-a.h.W- 


In  general,  Sq  >5^,  therefore,  Q^Qo,  and  eqn.  (26)  becomes 
^  =  0.5-a,^^ 

When  L  =  Ax.  h  =  .  and  qjc  =  ph.  eqn.  (27a)  becomes 


A'=  ^  1- 


CA-VSn 


(27a) 


(27b) 


Equation  (27b)  is  the  same  as  that  presented  by  Ponce  (1979)  and  Koussis 
(1976). 

MUSKINGUM  METHOD  WITH  VARIABLE  PARAMETERS 

Consider  an  {i+  l)th  period.  Let  /(/-)- 1)  be  the  inflow.  Qii+  1)  the  outflow. 
\\ii+  1)  the  storage,  and  K{i+  1)  the  corresponding  travel  time  for  that  period. 
Equation  (20)  can  be  rewritten  as 


In  the  same  way,  for  the  ith  period 
K(i)  = 

n  (I  +  2) 

By  dividing  eqn.  (28)  by  eqn.  (29).  the  following  is  obtained: 


m  = 


(29) 
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/<:(/+ 1) 

w(/-|- 1) 

”  QU) 

w(/+l) 

■  Q(0  T 

m 

w(i) 

.e(i+i)J 

w(0 

.Q{i+\)\ 

where  y  =  [a(2m+ 1 )]/(!+ a).  From  the  kinematic-wave  approximation 
(Ponce,  1979), 

Q  =  a.A^- 
Then 


e(/+l) 

ra2^(/-hl)Y^ 

rL^(/+i)i 

■w(/-f  i)Y^ 

Q{i) 

(X2A(i)  _ 

LAU)  J 

-  HO  - 

(31) 


By  substituting  eqn.  (31)  into  eqn.  (30) 

/«:(/+ 1)  _  W(/-M)r  VV(/)  y(2m+l)ft/(l+a) 

m  ~  w{i)  U(/4-l)J 

Three  interesting  cases  result  from  eqn.  (32): 

(1)  When  a  =  0.5,  m  =  1.5,  Pj  =  0.75,  [a(2m-}- l)/?2]/[l +«]  =  1-0,  eqn. 
(32)  simplifies  to 


HV+J) 

L  m;(/) 


(32) 


/*:(/+ 1) 

m 


=  1.0 


(33) 


Equation  (33)  expresses  that  the  travel  time  K  is  independent  of  storage.  In 
other  words,  vv(/-|- 1)  =  /r[A'(/-H  1 )/(/-!- 1) -I- [1 —A'(/-t- 1)]<2(^+ !)]■ 

(2)  When  [a(2m+ l)y52]/[l  4-a]<  1.0,  for  example,  a  =  0.5,  m  =  1.3, 
P.  =  0.75, 

p.a{2m+\)  ^ 

1+a 


KU+\) 

m 


rw(/+i) 


lO.IO 


L  w(i)  J 


(34) 


According  to  eqn.  (34),  the  travel  time  K  is  directly  proportional  to  storage  w. 
In  practice,  when  the  channel  cross-section  approximates  a  rectangular  shape 
and  the  channel  roughness  increases  with  the  increase  in  water  level,  the 
average  velocity  of  reach  travel  decreases  with  the  increase  in  the  water  level, 
and  the  travel  time  K  increases  with  the  increase  in  storage. 

(3)  When  [)52a(2w-|- !)]/(! -fa)  >  1.0,  for  instance,  a  =  0.6,  m  =  1.6, 
Pz  =  0.75, 

p.cc(2m+l) 


K(i+\) 

m 


n(/+l) 


-1-0.18 


L  nii)  J 


(35) 
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Equation  (35)  shows  that  travel  time  K  is  inversely  proportional  to  storage  u-. 
In  practice,  when  the  cross-section  approximates  a  triangular  or  parabolic 
shape,  channel  roughness  decreases  with  the  increase  of  the  water  level. 
From  eqn.  (27a),  the  following  formula  can  be  obtained: 

-''('■+1)  =  0.5-^[e(,+  l)]"  (27c) 

m  =  (27d) 

On  subtracting  eqn.  (27d)  from  eqn.  (27c) 

X(i+  1)  =  A-(/)+^{(e(,')f  _[2(/  +  ])](!>  pie) 

It  is.  thus,  seen  that  the  relationship  between  the  travel  time  and  storage 
depends  on  the  hydraulic  characteristics  of  the  channel. 

Three  versions  of  the  Muskingum  method  with  variable  parameters  are 
now  developed. 

Version  1:  dependence  of  travel  time  on  storage 


Equation  (32)  can  be  rewritten  as 


^(/+1)  =  K{i) 


L  w{i)  _ 


=  m 


H-(/+l) 

L  »»•(/■)  J 


(32a) 


where  y,  —  1.0  —  [a(2m  -I- 1  ))52]/(  1  +  x).  Taking  the  finite  difference  approxima¬ 
tion  of  eqn.  (8).  the  water  balance  equation  can  be  obtained  as 


u-(/+l) 


/(i+ I)-h/(i)  Q(i-h  l)  +  Q(i)  ^ 

- - At - r — ^-At  +  ir(/) 


(36) 


When  inflow  /(/+ 1),  u(/)  and  Q(/)  are  given.  0(/-t-l)  can  be  calculated  as 
follows: 

(1 )  Assume  a  value  of  ^"’(/-f  1). 

(2)  Calculate  »■(/+  1)  from  eqn.  (36)  using  2"’(/+  1). 

(3)  Calculate  AT/-I-1)  from  eqn.  (27e)  using  ^"’(/-i-l). 

(4)  Substitute  the  values  of  A:(/+1)  and  T(/-Hl)  into  eqns.  (2)-(5)  to 
calculate  the  Q(i+  1). 

(5)  Compare  Q(/+  1 )  with  Q‘(i+  I).  If  the  difference  between  Q(i-h  1 )  and 
0"'(/-i-l)  is  less  than  a  pre-assigned  error,  then  ^(i+l)  is  accepted  If 
(2(/+  1)>(2"'(/+  1).  then 
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Otherwise 


g">(/+l)-g(/4-l) 

2 


(6)  Substitute  2®(/+l)  for  return  to  step  (2)  and  continue  to 

iterate,  until  an  acceptable  value  of  Q(i+  1)  is  obtained. 


Version  2:  Dependence  of  travel  time  on  discharge  and  storage 


Ponce  (1979)  and  Koussis  (1976)  proposed  the  Muskingum  method  with 
variable  parameters.  The  travel  time  can  be  expressed  as 

^  A.x  ^x  ^Aa' 

"  C  ■  m)IA  ~  PQ 
With  Aa  =  L,  eqn.  (6a)  becomes 


A2sx  _  At  _  H’ 

iQ  ^JQ 


(6b) 


When  [a(2m+  !)]/(!  +a)  =  1.0,  eqn.  (20)  reduces  to  eqn.  (6b).  Therefore,  eqn. 
(6b)  is  a  special  case  of  eqn.  (20).  The  following  procedure  can  be  used  to 
obtain  the  routed  outflow: 

(1)  Assume  a  value  of  0*'*(/+l). 

(2)  Calculate  n’(/+l)  from  eqn.  (36). 

(3)  Estimate  A^(/+  1)  from  eqn.  (30)  using  ^•”(/+ 1)  and  vv(/+  1). 

The  rest  of  the  steps  are  the  same  as  for  the  first  method  described  above. 


Version  3:  dependence  of  travel  time  on  discharge  and  storage 


Equation  (30)  can  be  rewritten  as 

■]3t(2m+  1 )/( 1  +a) 

a:(/+i)  =  m\ . 


H’(/+  1) 

Qii)  7 

.  H’(/)  . 

(30) 


The  procedure  for  calculating  Qii+  1)  is  similar  to  that  in  the  first  case.  First, 
0"’(/+l)  is  assumed.  Then  *t’(/+l)  is  calculated  from  eqn.  (36),  Ar(/-l-l)  is 
estimated  from  eqn.  (30)  and  so  on. 


VALIDATION  OF  VARIABLE  PARAMETER  VERSIONS 


Inflow-outflow  data 


The  three  versions  of  the  Muskingum  method  with  variable  parameters 
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TABLE  1 


Parameters  of  the  three  versions  for  three  data  sets 


Version 

Parameters 

Data  set 

I 

II 

III 

1 

^(1) 

8.64 

33.19 

26.93 

A'd) 

0.13 

0.285 

0.285 

/T 

0.4 

-0.06 

-0.37 

0.001 

0.0014 

0.0005 

P 

0.5 

0.5 

0.5 

«■(!) 

1260 

450000 

450000 

2 

a:(1) 

16.15 

23.74 

18.97 

Y(l) 

0.05 

0.275 

0.275 

/ 

1.0 

1.0 

1.0 

3ti 

0.001 

0.0014 

0.0005 

P 

0.5 

0.5 

0.5 

M(l) 

1260 

450000 

450000 

3 

^:(i) 

10.88 

26.15 

22.51 

Y(l) 

0.05 

0.275 

0.275 

/ 

O.o 

1.15 

1.47 

0.001 

0.0014 

0.0005 

P 

0.5 

0.5 

0.5 

u(l) 

1260 

450000 

450000 

were  applied  to  three  actual  sets  of  inflow-outflow  data  (Table  1 ).  The  first 
data  set  was  taken  from  Wu  et  al.  (1985),  and  was  utilized  for  evaluating  the 
optimal  value  of  x  using  a  statistical  /-test  approach.  The  second  data  set  was 
taken  from  Wang  et  al.  (1987)  and  the  third  data  set  from  Yang  (1988).  Both 
the  second  and  third  data  sets  were  provided  by  the  Yangtze  River  Planning 
Office  in  China.  The  Yangtze  River  is  known  for  its  rampaging  ffoods  which 
have  caused  damage  of  catastrophic  proportions  throughout  history.  The 
river  reach  under  consideration  is  between  Wan  Xian,  Sichuan  Province,  at 
the  upstream  end  and  Yichang,  Hubei  Province,  at  the  downstream  end.  The 
three  data  sets  are  shown  in  Tables  2.  3  and  4,  respectively. 

Determination  oj  parameters 

The  three  versions  of  the  Muskingum  method,  discussed  above,  are  based 
on  the  nonlinear  storage  equation,  which  involves  five  common  parameters: 
u(l),  A'll).  ^(1).  a,  and  p.  In  addition,  the  first  version  also  includes  an 
exponent  •/,  =  \.Q  —  [y.{2m-\-\)(U]l{\-\-‘x):  and  the  third  version  includes 
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TABLE  2 


Muskingum  method  with  variable  parameters  for  data  of  Wu  et  al.  (1985) 


Recorded  Routed 


Time 

(h) 

Inflow 

(cfs) 

Outflow 

(cfs) 

m 

w(0 

(cfsh  ') 

wV) 

(cfsh  ') 

QU) 

(cfs) 

TS 

Q(l) 

(cfs) 

W 

QiO 

(cfs) 

TE 

0 

56 

70 

— 

1260.0 

_ 

_ 

.. 

6 

66 

66 

11.15 

1221.9 

723.1 

64.1 

66.6 

65.6 

12 

250 

102 

12.81 

1728.0 

1335.0 

82.6 

74.0 

48.0 

18 

550 

185 

16.82 

3418.2 

3430.5 

154.0 

135.0 

68.9 

24 

595 

265 

20.48 

5587.5 

6303.4 

267.9 

251.0 

183.1 

30 

420 

335 

22.12 

6773.0 

7963.5 

351.9 

338.4 

301.4 

36 

295 

370 

22.10 

6759.9 

7930.8 

367.5 

355.6 

342.6 

42 

210 

368 

21.26 

6136.9 

6995.0 

345.4 

335.4 

339.0 

48 

147 

310 

19.99 

5261.7 

5690.4 

303.6 

298.3 

313.3 

54 

100 

245 

18.50 

4334.2 

4327.6 

252.6 

254.7 

276.9 

60 

74 

200 

16.98 

3498.8 

3128.4 

199.9 

211.0 

235.8 

66 

60 

165 

15.63 

2841.1 

2212.5 

153.4 

172.8 

197.3 

72 

51 

132 

14.52 

2363.7 

1573.4 

116.8 

141.4 

164.4 

78 

46 

100 

13.67 

2035.2 

1150.0 

89.7 

116.3 

136.9 

TS;  for  the  results  obtained  from  this  study;  SM  =  2730.4,  A/  =  6  h. 

W;  for  the  results  obtained  from  Wu  et  al.  (1985);  SM  =  5525.6  where  X  =  0.1,  =  20.9  h, 

Ar  =  6h. 

TE;  for  the  results  obtained  by  graphical  procedure  of  trial  and  error;  SM  =  31  566.0  where 
X  =  0.2.  K  =  26.7  h,  Ar  =  6h. 

y  =  [a(2m+  1)]/(1  +a).  If  the  hydraulic  characteristics  of  the  channel  reach 
are  known,  then  the  exponent  yi  or  y  can  be  estimated  from  a,  and  m; 
a,  =  (boli)lilLSQ)  and  ^  can  be  determined  from  the  stage-discharge  curve  of 
steady  flow  of  the  channel.  The  initial  storage  of  reach  w(l)  can  be  determined 
from  the  initial  inflow,  outflow  and  the  length  or  the  corresponding  travel  time 
of  the  reach.  Therefore,  the  parameters  to  be  determined  in  the  three  versions 
are  similar  to  those  in  the  conventional  Muskingum  method,  that  is,  Af(l)  and 
A"(l).  If  the  hydraulic  characteristics  of  the  channel  reach  are  unknown,  then 
y, ,  or  y  are  to  be  estimated  from  observed  inflow  and  outflow  hydrographs. 
Since  the  parameters  to  be  estimated  are  only  a  few,  the  direct  search  method 
was  used  to  obtain  their  optimum  values  by  minimizing  the  sum  of  the  squares 
of  differences  between  observed  and  computed  outflows. 

The  direct  search  methods  require  evaluation  of  only  the  objective  function 
and  do  not  use  the  partial  derivatives  of  the  function  in  finding  the  minimum. 
These  methods  are  most  suitable  for  simple  problems  involving  a  relatively 
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TABLE  3 


Muskingum  method  with  variable  parameters  for  data  of  Wang  et  al.  (1987) 


Recorded 

Routed 

Time 

Inflow 

Outflow 

X{t) 

m 

w(t) 

w'{t) 

Q(t) 

Q(t) 

Q(t) 

(m^s”') 

(m's-') 

(h) 

(m^s-'h- 

')  (m's-'h- 

')  (m's-') 

(m's-‘) 

(m's-‘) 

(10') 

(10') 

SV 

W 

TE 

7.27 

17900.0 

17900.0 

0.24 

23.74 

454.8 

424.1 

17900.0 

17900.0 

17900.0 

28 

18200.0 

18000.0 

0.239 

26.13 

609.0 

567.9 

18000.0 

18000.0 

17994.8 

29 

25700.0 

18100.0 

0.233 

29.17 

879.37 

812.9 

18952.0 

19287.0 

20494.9 

30 

36500.0 

23700.0 

0.205 

30.11 

1008.3 

919.4 

24878.6 

25079.0 

27198.5 

31 

38700.0 

32700.0 

0.172 

27.39 

993.7 

904.6 

33264.8 

32521.0 

33773.4 

8.1 

35700.0 

36300.0 

0.161 

26.17 

899.7 

818.9 

36388.4 

35641.0 

35939.4 

2 

30500.0 

35400.0 

0.168 

26.29 

841.6 

764.8 

34493.7 

34486.0 

34144.6 

3 

28400.0 

31800.0 

0.182 

27.28 

962.86 

874.8 

31000.1 

31358.0 

31177.1 

4 

35300.0 

30200.0 

0.184 

29.31 

1080.8 

979.5 

30037.6 

30510.0 

31602.4 

5 

39100.0 

35100.0 

0.168 

28.87 

1125.3 

1018.1 

34007.9 

33989.0 

35140.7 

6 

39200.0 

38300.0 

0.157 

28.05 

1130.76 

1022.96 

37299.5 

36963.0 

37675.0 

7 

38300.0 

38400.0 

0.154 

28.00 

1111.44 

1005.5 

38089.6 

37911.0 

38354.5 

8 

36800.0 

37000.0 

0.156 

28.17 

1096.7 

992.07 

37665.4 

37580.0 

37846.0 

9 

36000.0 

36300.0 

0.159 

28.59 

1337.2 

1209.7 

3667.8 

36644.0 

36932.0 

10 

48300.0 

38800.0 

0.153 

30.10 

1770.3 

1589.2 

38135.7 

38117.0 

40230.3 

11 

65400.0 

47900.0 

0.123 

29.39 

1993.3 

1767.4 

47886.9 

47205.0 

50735.8 

12 

67700.0 

58400.0 

0.0896 

27.39 

1888.4 

1672.4 

59521.0 

58890.0 

60732.1 

13 

56800.0 

61000.0 

0.0811 

26.65 

1692.9 

1497.5 

61968.9 

62204.0 

61691.6 

14 

47500.0 

56300.0 

0.097 

27.34 

1457.5 

1284.5 

56045.3 

56789.0 

55660.2 

15 

37800.0 

47800.0 

0.122 

28.48 

1227.6 

1075.4 

47986.7 

48646.0 

47436.3 

16 

29600.0 

38900.0 

0.150 

29.86 

1043.4 

907.4 

39616.0 

39918.0 

38753.4 

17 

23900.0 

31500.0 

0.179 

31.68 

964.5 

834.4 

32001.7 

32098.0 

31165.9 

18 

23200.0 

26000.0 

0.199 

33.91 

1147.5 

992.2 

26696.6 

26612.0 

26351.8 

19 

33100.0 

24500.0 

0.204 

38.50 

1491.4 

1281.3 

25116.7 

25954.0 

27488.5 

20 

45200.0 

30300.0 

0.181 

39.81 

1736.98 

1472.2 

30262.5 

32445.0 

34859.9 

21 

49500.0 

39500.0 

0.150 

36.89 

1715.7 

1448.5 

39070.6 

40979.0 

42755.0 

22 

42900.0 

43800.0 

0.134 

34.56 

1583.9 

1337.0 

43938.8 

44802.0 

44932.6 

23 

36300.0 

41200.0 

0.141 

34.50 

1514.8 

1276.5 

42119.8 

42027.0 

41561.8 

24 

35300.0 

38600.0 

0.154 

35.78 

1551.2 

1306.7 

38396.8 

37883.0 

37919.7 

25 

38200.0 

36100.0 

0.159 

37.25 

1518.7 

1279.3 

36668.1 

36463.0 

37180.3 

26 

31300.0 

34900.0 

0.165 

36.95 

1313.4 

1104.6 

35466.7 

34946.0 

34713.9 

28 

27500.0 

31800.0 

0.176 

37.72 

1219.0 

1023.3 

32443.9 

31724.0 

31354.9 

29 

25000.0 

28700.0 

0.189 

38.88 

1183.0 

991.64 

29237.9 

28419.0 

28128.0 

SV:  the  results  obtained  from  the  special  version  yield  SM  =  15  507000.0,  At  =  24  h. 

W:  the  results  obtained  by  Wu  et  al.  (1985),  SM  =  26  21 1  660.0.  X  =  QM.  K  =  30.9  h. 
TE;  the  results  obtained  by  graphical  method  of  trial  and  error  yield  SM  =  88  740000  0 
A' =  0.0. /r  =  26.0  h. 
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TABLE  4 


Muskingum  method  with  variable  pararameters  for  data  of  Yang  et  al.  (1988) 


Recorded 

Routed 

Time 

Inflow 

Outflow 

TV 

0 

TE 

m 

Qit) 

7.2 

21300.0 

22200.0 

22.67 

22000.0 

22000.0 

22200.0 

0.14 

27900.0 

22800.0 

27.78 

21646.9 

22680.0 

22962.0 

3.2 

38800.0 

23400.0 

24.90 

23671.4 

26700.0 

27693.0 

14 

47300.0 

30100.0 

20.63 

30680.3 

33900.0 

35305.0 

4.2 

52000.0 

39900.0 

18.03 

39825.3 

41170.0 

42678.0 

14 

53800.0 

46800.0 

16.84 

47049.2 

46760.0 

48057.0 

5.2 

52700.0 

50900.0 

16.64 

50858.6 

50080.0 

50953.0 

14 

48900.0 

51100.0 

17.05 

51191.6 

50830.0 

51179.0 

6.2 

43400.0 

49200.0 

18.10 

48930.6 

49110.0 

48906.0 

14 

37800.0 

44800.0 

19.36 

44751.8 

45300.0 

44871.0 

7.2 

32000.0 

39900.0 

20.86 

40127.3 

40890.0 

39951.0 

14 

26900.0 

34900.0 

22.77 

35298.6 

35790.0 

34688.0 

8.2 

23400.0 

30600.0 

24.64 

30638.4 

30920.0 

29817.0 

14 

21400.0 

26700.0 

26.21 

26923.2 

26980.0 

25970.0 

9.2 

19600.0 

23700.0 

27.67 

24168.4 

23950.0 

23159.0 

14 

18600.0 

21800.0 

28.02 

22000.0 

21670.0 

21046.0 

10.2 

17200.0 

20100.0 

28.69 

20675.4 

19960.0 

19459.0 

14 

16400.0 

18900.0 

29.64 

19288.2 

18490.0 

18087.0 

11.2 

16400.0 

17800.0 

30.40 

18101.2 

17460.0 

17175.0 

14 

17000.0 

17300.0 

17434.0 

17020.0 

16870.0 

TV:  the  results  from  this  version,  SM  =  3  004235.0,  At  =  12h. 

O:  the  results  obtained  from  the  method  of  optimizing  X  value  yield  SM  =  31  225  100.0, 
X  =  0.123,  K  =  20.92  h.  At  =  12  h. 

TE:  the  results  obtained  by  graphical  method  of  trial  and  error  yield  SM  =  58656232.0, 
A'  =  0.10,  K  =  18  h. 


small  number  of  parameters.  All  the  direct  search  methods  are  iterative  in 
nature  and  hence  they  require  an  initial  point  Yi  to  start  the  iterative 
procedure,  and  differ  from  one  another  only  in  the  method  of  generating  the 
new  point  (from  Y^)  and  in  testing  the  point  7,+  ,  for  optimality.  In  this 
study,  the  simplex  method  was  used  to  optimize  the  parameters.  The  basic 
idea  of  this  method  is  to  compare  the  values  of  the  objective  function  at  the 
n  + 1  vertices  of  a  general  simplex  and  move  this  simplex  gradually  towards 
the  optimum  point  during  the  iterative  process.  The  movement  of  the  simplex 
is  achieved  by  using  three  operations  known  as  reflection,  contraction  and 
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expansion.  A  complete  discussion  of  the  simplex  method  can  be  found  in 
Himmelblau  (1972). 

According  to  the  general  principle  of  determining  the  parameters 
mentioned  above,  the  three  versions  of  the  Muskingum  method  were  applied 
to  the  three  data  sets.  The  initial  storage  h'(I)  was  estimated  from  the  initial 
inflow,  outflow  and  the  reach  length  or  the  average  travel  time  of  the  reach. 
The  rating  curve  for  the  steady  flow  was  used  to  determine  jS  and 
a,  =  (^?o/?)/(2L5'o).  Then  the  direct  search  method  was  used  to  optimize  A^(l), 
X(l),  7 1  or  y.  For  the  second  (special)  version,  y  =  1.0,  only  /r(l)  and  A'(l) 
needed  to  be  optimized.  The  parameters  of  the  three  versions  for  three  data 
sets  are  listed  in  Table  1  which  shows  that  the  initial  storage  vv(l)  is  the  same 
in  the  three  versions  developed.  The  values  of  ^'(l),  estimated  using  a  direct 
search  method,  are  the  same  except  for  the  first  version. 

Flow  routing 

Each  version  with  the  parameters  obtained  above,  as  well  as  the  graphical 
trial  and  error  method,  were  applied  to  all  three  data  sets.  For  purposes  of 
comparison  of  the  methods,  the  sum  {SM)  of  the  squares  of  the  deviations  {d) 
between  observed  and  computed  outflows,  SM  =  was  utilized. 

Application  of  version  1 

The  results  of  flow  routing  by  the  first  version  are  shown  only  for  the  first 
example  data  set  to  save  space.  Table  2  includes  recorded  inflow  and  outflow, 
as  well  as  routed  outflow  obtained  from  the  first  version  and  the  graphical 
method.  Also  included  are  the  results  of  Wu  et  al.  (1985).  The  results  of 
calculations  show  that  for  At  =  6h,  SM  =  2730.4  for  this  version. 
SM  =  5523.6  for  the  method  of  Wu  et  al.  (1985)  with  A'  =  0.1  and 
K  =  20.9  h  and  SM  =  31  566.0  by  the  graphical  method  with  X  =  0.2  and 
K  =  26.7  h.  Clearly  for  this  data  set  the  first  version  of  the  Muskingum 
method  was  the  best  of  the  three  methods. 

Application  of  version  2  (special  case) 

The  results  of  flow  routing  by  the  special  version  are  shown  only  for  the 
second  example  data  set.  Table  3  includes  inflow  and  outflow,  as  well  as  the 
routed  outflow,  3^(1)  and  A!'(l)  obtained  from  the  special  version  and  the 
graphical  method.  Also  included  are  the  results  obtained  from  the  method  of 
Wu  et  al.  (1985).  The  computed  results  show  that  for  A/  =  24  h. 
SM  =  15  507  000.0  for  this  version,  SM  =  26211  660.0  for  the  method  of 
Wu  et  al.  (1985)  with  A'  =  0.170  and  =  30.9  h  and  5'M  =  88  740000.0  by 
the  graphical  method  with  A'  =  0.0,  A:  =  26.0  h.  For  this  data  set.  the  special 
version  of  the  Muskingum  method  was  the  best  of  the  three  methods. 
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TABLE  5 


Comparison  of  the  results  obtained  from  the  three  versions  developed  with  those  of  the 
conventional  trial-and-error  method 


Event 

Routed  outflow  (fd') 

First  version 

Second  version 
(special) 

Third  version 

Conventional 

trial-and-error 

method 

1 

2730.4 

8.7% 

4714.8 

14.9% 

3866.0 

12.2% 

31566.0 

100% 

2 

17101110.0 

19.25% 

15507000.0 

17.5% 

13811000.0 

15.56% 

88740000.0 

100% 

3 

8484371.0 

14.46% 

9119755.0 

15.55% 

3004235.0 

5.12% 

58656232.0 

100% 

Application  of  version  3 

The  results  of  calculation  for  the  third  version  are  shown  only  for  third  data 
set  of  Yang  (1988).  The  routed  results  are  given  in  Table  4.  The  flow  was  also 
routed  using  the  method  of  optimizing  X  as  well  as  the  graphical  method. 
SM  =  3  004235.0  for  this  version,  =  31  225  100.0  for  the  method  by  Yang 
(1988)  with  X  =  0.123,  K  =  20.92  h  and  SM  =  58656232.0  by  the  graphical 
method  with  Y  =  0.1,  Y  =  18h,  A/  =  1 2  h.  Clearly,  for  this  data  set,  the  third 
version  of  the  Muskingum  method  was  the  most  accurate  of  the  three  methods. 

Discussion  of  the  results 

The  three  versions  of  the  Muskingum  method  with  variable  parameters 
were  compared  with  the  graphical  method  for  the  three  data  sets.  The  method 
presented  by  Ponce  (1979)  and  Koussis  (1976)  is  regarded  as  a  special  case  of 
the  method  developed  in  this  study.  The  results  of  this  version  were  also 
compared  with  the  other  two  versions,  and  are  summarized  in  Table  5  which 
shows  that  the  methods  developed  in  this  study  are  more  accurate  than  the 
conventional  trial  and  error  method.  For  the  third  event,  the  squared  error 
obtained  from  the  third  version  was  only  5.12%  of  that  of  the  conventional 
trial  and  error  method.  Even  for  the  worst  of  the  three  versions  for  three 
events,  the  error  was  19.25%  of  that  of  the  conventional  trial  and  error 
method.  Table  5  also  shows  a  comparison  of  the  results  from  the  special 
version  (second)  presented  by  Ponce  (1979)  and  Koussis  (1976)  with  those 
from  the  other  two  versions.  For  the  first  data  set.  the  first  version  gave  the 
best  result;  for  the  second  and  third  data  sets,  the  third  version  produced  the 
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best  results.  Only  for  the  second  data  set  was  the  result  from  the  second 
version  better  than  the  first  version. 

It  is  well  known  that  the  routing  accuracy  of  the  conventional  Muskingum 
method  depends  on  the  relationship  between  the  storage  w{i)  and 
Q'iO  [Q'{^  =  +  If  an  optimum  method  is  used  to  obtain  the 

value  X  which  makes  the  relationship  between  the  storage  w{i)  and  Q'{i)  a 
straight  line,  then  the  routing  result  is  the  best  that  can  be  achieved  by  the 
conventional  method;  this  is  possible  only  under  the  condition  that 
[a(2m+l)i?2]/(l+a)  =  1.0. 

In  general,  [a(2/M+  \)^2]l{\  +a)  is  not  equal  to  1.0,  so  that  it  is  impossible 
to  optimize  the  value  of  X  to  make  the  relation  between  w{i)  and  Q'{i)  a 
straight  line.  In  some  channel  reaches,  the  relation  w(/)  and  Q'{i)  is  a  curve 
in  which  case  the  reach  travel  time  K  is  not  constant. 

The  parameters  K  and  X  vary  with  time  in  this  study.  It  is  possible  to  make 
the  relationship  between  the  storage  calculated  from  the  water  balance 
equation  and  that  calculated  from  the  storage  equation  h’'(/-|-1)  = 
Ar(/+ 1){A'(/+ 1)/(/-|- 1)  +  [I  —  XCi-f  1)]^(/+ 1)}  plot  as  a  straight  line. 

Figure  1  shows  the  relation  between  w  and  vv'  for  the  first  flood  event, 
where  the  ordinate  is 


=  /('  + 1 ) + A-)  _  g(- + 1 ) + 0(0  ^ 


and  the  abscissa  is  1)  =  K(i+  1){A'(/+  l)/(/+ 1)  +  [1  -X(i+  \)]Q{i+ 1)} 
where 


/:(/+!)  =  m 


1) 

.  MO 


It  can  be  seen  that  the  data  points  from  the  first  version  are  quite  close  to  a 
straight  line.  The  results  from  the  conventional  Muskingum  method  are  also 
shown  in  Fig.  1,  where  the  loop  is  formed  by  the  data  points. 

Figure  2  shows  the  relation  between  w  and  w'  for  the  second  flood  event 
of  Changjiang  River  (Yangtze),  where  the  abscissa  is 


h'(/+I)  =  /!:(/+ l){Y(/+l)/(/+l)  +  [l-Y(/+l)]2(/+l)} 
where 


A:(/+1)  =  m 


w(/+  !)■ 

Q(0  1 

_  MO  - 

Le(/+i)J 

It  can  be  seen  that  the  data  points  almost  lie  on  the  straight  line  (if  all  the 
points  lie  exactly  on  a  straight  line,  then  the  computed  outflows  will  be  just 
the  same  as  the  observed  outflows).  The  relation  of  w  with  u  '  calculated  from 
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Fig.  1.  Relationship  between  channel  storage  and  m'  =  F{K,X,l,Q)  for  event  I. 


the  conventional  Muskingum  method  is  also  shown  in  Fig.  2,  where  a  loop  is 
formed  by  the  points. 

Figure  3  shows  the  third  flood  event  of  Changjiang  River,  where  the 
abscissa  is  vvX/+l)  =  ^(/+ 1){X(/+ !)/(/+ 1)  +  [1 -A'(/+ 1)]0(/+ 1)}  where 


K{i+\)  =  m 


w(i+  1) 

(2(0  1 

-  W'(/)  _ 

L0(/+i)j 

It  can  be  seen  that  the  data  points  almost  lie  on  the  curve.  The  results  from 
graphical  trial  and  error  methods  are  also  shown  in  Fig.  3,  where  the  shape 
of  the  curves  formed  by  the  data  points  calculated  from  the  two  methods  are 
similar,  in  which  case  the  reach  travel  time  is  not  constant;  therefore,  the  discharge 
of  outflows  estimated  from  this  study  are  more  accurate  than  the  conventional 
Muskingum  method  in  which  the  reach  travel  time  is  considered  constant. 


SUMMARY  AND  CONCLUSIONS 


In  this  study,  three  versions  of  the  Muskingum  method  with  variable 
parameters  were  derived  by  simplifying  the  St.  Venant  equations  for  unsteady 
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W  =  F  (K.X,  !.Q)  (C«S  OOOO*  HRS) 


Fig.  2.  Relationship  between  channel  storage  and  w  =  F(K.XJ,Q)  for  event  2. 


flow.  These  versions  are  simple  to  apply,  and  it  is  relatively  easy  to  determine 
the  parameters  using  a  direct-search  method.  The  following  conclusions  are 
drawn  from  this  study. 

(1)  The  relationship  between  the  reach  travel  time  and  storage  depends  on 
the  hydraulic  characteristics  of  the  channel  reach.  If  [a(2m-l-  l)j52]/(l  -l-a)<  1.0, 
then  the  relation  between  the  travel  time  and  storage  is  in  direct  proportion; 
if  [a(2m-f  l)iS2]/(l  -f  a)>  1.0,  then  it  is  in  inverse  proportion;  if  [a(2m-|- 1)/?2]/ 
( 1  +  a)  =  1.0,  a  special  case  of  this  relation  does  not  exist.  These  relations 
between  storage  and  travel  time  are  identified  by  routing  observed  floods. 

.(2)  The  Muskingum  method  with  variable  parameters  can  overcome  the 
drawback  of  the  conventional  graphical  trial  and  error  loop  method  or  of 
optimum  methods,  which  are  inadequate  to  represent  the  storage  hysteresis. 

(3)  The  three  versions  of  the  Muskingum  method  are  several  times  more 
accurate  than  the  conventional  trial  and  error  method. 

(4)  The  reach  travel  time  which  depends  on  the  storage,  channel  charac¬ 
teristics  and  discharge  is  obtained  from  the  simplification  of  the  St.  Venant 
equations,  when  the  exponent  of  discharge  [ci(2m+  1)^2 ]/(l  +«)  =  1.0,  eqn. 
(20)  reduces  to  eqn.  (6)  presented  by  Ponce  (1979)  and  Koussis  (1976).  which 
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Fig.  3.  Relationship  between  channel  storage  and  w  =  F(K,XJ,Q)  for  event  3. 


is  regarded  as  a  special  case  of  the  method  developed.  The  results  show  that 
the  special  version  is  several  times  more  accurate  than  the  conventional  trial 
and  error  method,  but  less  accurate  than  the  other  two  versions. 
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ABSTRACT 

Wang,  G.-T..  Singh,  V.P.  and  Yu,  F.X.,  1992.  A  rainfall-runoff  model  for  small  watersheds.  J.  Hydro!., 
138:  97-117. 

A  rainfall-runoff  model  was  developed  by  combining  the  excess-rainfall  process  and  the  runoff- 
concentration  process.  The  excess  rainfall  was  modeled  by  using  the  two-parameter  Green-Ampt  infiltra¬ 
tion  approach.  A  six-parameter  linear-discrete  model  was  used  to  model  the  runoff  hydrograph.  The 
infiltration  parameters  were  estimated  by  using  the  simplex  method,  and  the  runoff  parameters  by  least 
squares.  The  model  was  calibrated  on  ten  watersheds  and  verified  on  seven.  The  model-simulated  runoff 
hydrographs  were  in  close  agreement  with  observed  runoff  hydrographs. 


INTRODUCTION 

For  the  purposes  of  modeling,  the  rainfall-runoff  process  is  divided  into  the 
excess-rainfall  process  and  the  runoff-concentration  process.  The  former 
focuses  on  estimating  the  infiltration  process,  and  the  latter  on  calculating  the 
outlet  hydrograph  from  excess  rainfall  using  a  mathematical  model. 

In  1911.  Green  and  Ampt  developed  an  infiltration  equation,  with  the 
horizontal  wetting  front  separating  the  wetted  and  dry  parts  of  the  soil  during 
soil-water  movement.  Many  infiltration  formulae  have  since  been  proposed 
(see  Singh  (1989)  for  a  recent  review).  The  Green-Ampt  (GA)  equation, 
however.' has  lately  been  receiving  a  great  deal  of  attention.  Mein  and  Larson 
(1973)  and  Swartzendruber  (1974)  developed  modified  versions  of  the  GA 
equation  for  determining  the  ponding  time  and  modeled  infiltration  during  a 
steady  rain.  Chu  (1978)  extended  their  analysis  to  unsteady  rainfall  events.  He 
used  two  time  parameters  (the  ponding  time  and  the  pseudotime)  to  modify 
the  GA  equation,  to  describe  infiltration  during  an  unsteady  rain. 

After  excess  rainfall  is  calculated,  the  runoff  hydrograph  can  then  be 
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estimated  using  a  suitable  runoff  model.  In  1957,  Nash  reported  a  two- 
parameter  linear  model  with  N  reservoirs  in  series.  The  Nash  model  has  been 
extensively  used  in  applied  hydrology.  Most  runoff  models  are  continuous¬ 
time  models  (Singh,  1988),  although  discrete-time  models  have  also  been 
developed.  O’Connor  (1982)  used  the  transfer  function  approach  to  derive 
discretely  coincident  forms  of  several  linear,  time-invariant  models.  Wang  and 
Yu  (1986)  suggested  a  more  general  linear  discrete  model  from  which  some  of 
the  existing  models  can  be  derived  as  special  cases.  These  models  attempt  to 
establish  a  relationship  between  excess-rainfall,  calculated  using  an  infil¬ 
tration  formula,  and  the  resulting  direct  runoff  hydrogaph. 

Mays  and  Taur  (1982)  used  a  nonlinear  programming  formulation  to 
determine  simultaneously  the  coordinates  and  the  period  of  excess  rainfall. 
Wang  and  Yu  (1990),  and  Wang  et  al.  (1992)  modeled  the  rainfall-runoff 
process  by  considering  losses.  In  these  models,  the  number  of  interval  losses 
computed  was  equal  to  the  number  of  periods  of  rainfall  intervals  used.  In 
practical  applications,  however,  it  is  difficult  to  use  these  models  for  predicting 
the  runoff  hydrograph,  because  the  number  of  rainfall  excess  intervals  of  one 
event  may  not  be  exactly  equal  to  that  for  another.  Therefore,  the  computed 
interval  losses  cannot  be  used  for  predicting  another  event. 

In  this  study,  an  effort  was  made  to  combine  the  excess-rainfall  process  and 
the  runoff-concentration  process  into  a  unified  process.  The  computation  of 
excess  rainfall  during  an  unsteady  rain  has  been  given  by  Chu  (1978).  With  the 
use  of  measured  rainfall  hyetographs  and  runoff  hydrographs,  optimal  values 
of  infiltration  and  runoff  model  parameters  were  simultaneously  estimated. 

RAINFALL-RUNOFF  PROCESS  DURING  UNSTEADY  RAIN 


Determination  of  excess  rainfall 


For  many  rainfall  events,  there  is  an  initial  period  during  which  all  of  the 
rainfall  infiltrates  into  the  soil.  During  this  period,  the  capacity  of  the  soil  to 
absorb  water  decreases  until  it  becomes  less  than  the  rainfall  intensity.  At  this 
point,  the  ground  surface  becomes  ponded  with  water.  As  rainfall  continues, 
the  surface  pondage  exceeds  the  surface  retention  capacity  and  the  runoff 
begins.  Under  ponded  conditions,  the  infiltration  process  is  independent  of 
the  time  distribution  of  rainfall.  For  an  unsteady  rainfall  event,  there  may  be 
several  periods  during  which  the  rainfall  intensity  exceeds  the  current  infil¬ 
tration  rate  and  ponding  may  appear  and  disappear. 

The  GA  equation  describes  the  infiltration  process  under  a  ponded  surface, 
and  can  be  written  (Mein  and  Larson.  1973)  as 


1  + 


F 


(1) 
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where /p  is  the  infiltration  capacity  (mm min"'),  K  is  the  average  hydraulic 
conductivity  of  the  wetted  zone  (mm  min“' ),  S  (mm)  and  M  are  the  difference 
of  the  average  capillary  potential  and  the  difference  of  the  average  soil 
moisture  before  and  after  wetting,  respectively,  expressed  as  volume  fraction 
of  soil  (volume/volume),  and  F  is  the  cumulative  infiltration  (mm). 

If  we  let  /p  be  the  time  at  which  surface  ponding  starts  from  a  stage  without 
surface  ponding,  be  the  cumulative  infiltration  at  /  =  tp,  and  F^  be  the 
cumulative  infiltration  under  ponding,  integration  of  (1)  from  t  =  t^io  t  and 
from  F  =  Fq  to  Fp,  and  rearrangement,  yield. 


A 

SM 


In  1  + 


SM 


+  In 


K{t  -  tp) 

SM 


(2) 


In  practice,  the  duration  of  a  rainfall  event  is  divided  into  many  short  periods 
in  such  a  way  that  within  each  period,  the  rainfall  intensity  is  essentially 
constant.  For  such  a  case. 


/(/)  =  =  /  =  constant,  «  =  0,  1,  2,  .  .  .  (3) 

tn-  tn-\ 


where  i(t)  is  the  rainfall  intensity,  t„  —  is  the  nth  time  interval  (short), 
F(/„)  is  the  cumulative  rainfall  in  millimeters  at  the  end  of  the  nth  time 
interval,  and  F(/„_i)  is  the  cumulative  rainfall  in  millimeters  at  the  beginning 
of  the  nth  time  interval. 

The  variable  cumulative  rainfall  P{t)  within  a  short  period  can  be  written 
as 


Pit)  =  F(/„_,)  +  I  /(r)d/  =  F(t„_,)  +  (/  -  t„_,)/  (4) 

Combining  (2),  (3)  and  (4),  and  letting  t  =  tp,  we  obtain  the  ponding  time  for 
this  special  case: 

fit)  =  /•(/)  <  /p  (5) 


dF 

d! 


0 


(6) 


Assuming  that  the  evaporation  during  a  rain  period  is  negligible,  the  infil¬ 
tration  process  during  an  unsteady  rain  can  now  be  separated  into  two  stages. 
In  the  first  stage,  there  is  no  surface  ponding  for  the  period  from  to  t„.  In 
this  period,  the  infiltration  rate  and  cumulative  rainfall  excess  F(0  satisfy 

R(t)  =  F(t„_,).  ^  ^  (7) 
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The  cumulative  infiltration  can  be  obtained  by  the  mass-balance  principle  as 

F{t)  =  P(t)  -  (8) 

In  the  second  stage,  there  is  surface  ponding  for  the  period  from  to  t„.  In 
this  period,  the  infiltration  rate  is  given  by  (1),  and  the  cumulative  rainfall 
excess  R(t)  can  be  expressed  using  the  mass-balance  principles; 


R{t)  =  P{t)  -  Fp  -  A  <  r  ^  7;  (9) 

where  D  is  the  surface  retention  capacity.  Let  us  define  a  time  constant  called 
pseudotime 


K 


t  — 


+  h 


SM 


(10) 

(11) 


Surface  ponding  occurs  when  rainfall  intensity  equals  the  infiltration 
capacity,  which  is  defined  as  the  rate  of  infiltration  that  reaches  its  maximum 
capacity  for  a  given  type  of  soil  and  moisture  condition,  so  that 

=  /p  (12) 

with  use  of  (12),  (1)  becomes 


,  _  KSM 
”  ■  /(/p)  -  K 


(13) 


Application  of  (8)  at  the  time  of  ponding  to  (13)  yields 


PU,)  -  Pin 


KSM 
i{t,)  -  K' 


i  >  K 


(14) 


Equation  (14)  is  an  implicit  expression  for  determining  the  ponding  time 

After  /p  is  determined,  the  pseudotime  can  then  be  calculated  by  (1). 
where 


Po  =  PUp)  -  R(n,  i„-i  <  t'  <  /„ 

For  a  constant  rainfall  intensity,  we  can  write 


at)  = 


P{t„)  -  P{t„_,) 


-  t. 


=  /  =  constant 


n  —  I 


(15) 

(16) 


The  ponding  time  can  be  obtained  simply  by  combining  (4)  and  (14)  and 
letting  t'  =  and  t  = 


'p  = 


KSM 
1  -  K 


—  P{t„-\)  +  RUn- I  ) 


+  In 


I 


(17) 
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Discrete  linear  model 

The  discrete  excess  rainfall-runoff  model  is  applied  to  calculate  the 
hydrograph  of  a  watershed  from  the  excess  rainfall  estimated  as  above.  The 
model  can  be  written  as 

Q{t)  =  a\Q{t  —  1)  +  —  2)  +  ...  +  OpQit  —  p) 

+  b^I{t)  +  b,I{t  -  1)  +  .  .  .  +  b^I{t  -  q)  (18) 

where  Q{t)  and  /(/)  are  the  direct  runoff  and  excess  rainfall  at  time  /,  respect¬ 
ively,  and  as  and  bs  are  the  time-invariant  parameters.  A  complete  discussion 
of  the  discrete  linear  model  has  been  given  by  Wang  and  Yu  (1986)  and  Wang 
et  al.  (1992). 

DETERMINATION  OF  MODEL  PARAMETERS 

The  rainfall-runoff  process  model  includes  two  groups  of  parameters.  In 
one  group  are  soil  infiltration  process  parameters  K,  S  and  M.  Parameters  S 
and  M  always  appear  together;  thus,  the  product  SM  can  be  treated  as  one 
parameter  which  represents  the  antecedent  soil  moisture  condition.  This 
reduces  to  two  the  number  of  parameters  for  the  soil  infiltration  process.  The 
number  of  parameters  for  excess  rainfall-runoff  process  depends  on  the  order 
of  the  model.  In  general,  six  parameters  of  excess  rainfall-runoff  process 
model  suffice  to  calculate  the  hydrograph.  Therefore,  for  the  rainfall-runoff 
process  model,  eight  parameters  are  to  be  estimated  from  measured  rainfall 
hyetographs  and  runoff  hydrographs. 

In  this  paper,  a  direct  search  method  was  used  to  estimate  the  infiltration 
parameters  (K  and  SM).  The  excess  rainfall  hyetograph  was  then  calculated, 
and  the  excess  rainfall-runoff  model  parameters  were  determined  by 
minimizing  the  sum  of  squares  of  differences  between  the  observed  and 
calculated  discharges.  The  parameter  estimation  method  can  be  described  in 
the  following  steps: 

(1)  A  rainfall  hyetograph  is  divided  into  several  periods  during  which  the 
rainfall  intensity  is  constant  and  the  value  of  retention  capacity  D  is  selected. 

(2)  Starting  from  the  first  period,  the  rainfall  intensity  /(/))  is  compared 
with  K.  If  /(/)  is  less  than  K.  then  excess  rainfall  R{t)  =  0  and  the  cumulative 
infiltration  F(t)  =  Pit).  If  /(/)  is  greater  than  K.  then  the  surface  condition 
indicator  can  be  estimated  as 

,  KSM  ,,Q, 

c.  =  PC)  -  RC  -  1)  -  05) 

If  Q  is  less  than  zero,  excess  rainfall /?(0  =  0,  andF(/)  =  P(/).  If  Q  is  greater 
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than  zero,  then  the  ponding  time  can  be  calculated  from  (17),  and  the 
cumulative  rainfall  prior  to  ponding  time  can  be  computed  as 


K 


tv  = 


SK 

I  -  K 


+  T„., 


(20) 


/>(/p)  =  />(/  -  1)  +  [/p  -  (/  -  l)A/]/(0  (21) 

and  we  shift  the  time  as 


P(tp)  -  /?(/  -  1)  ,  P(tp)  -  P(/  -  1)1  SM 

- ? - ‘"L‘  - SM - 


The  time  on  the  shifted  time-scale  is 


^  ^  ~  +  4 

and  Fp  can  be  calculated  as 


(22) 

(23) 

(24) 


This  formula  is  implicit  in  Fp,  so  the  Newton-Raphson  iteration  can  be  used 
to  estimate  Fp .  The  cumulative  excess  rainfall  is  calculated  as 


m  =  P(t)  -  Fit)  -  D 


(25) 


The  excess  rainfall  is  calculated  for  the  next  period,  and  so  on. 

(3)  After  the  excess  rainfall  is  estimated,  the  least-squares  method  can  be 
used  to  determine  the  parameters  of  the  excess  rainfall-runolf  model  as 


where 


0 

0 

0 

/(I) 

0 

0 

Q(\) 

0 

0 

1(2) 

/(I) 

0 

(2(2) 

01) 

0 

7(3) 

7(2) 

/(I) 

Q(5) 

(2(2) 

01) 

7(4) 

7(3) 

7(2) 

• 

• 

I(k) 

I(k  - 

1)  Kk  -  2) 

0 

0 

0 

Qim  - 

1)  Qim  - 

2)  Qim  -  3) 

(26) 


(27) 


E  =  [Qi\)Qi2)QO)  .  .  .  Qim)Y 


(28) 
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(4)  We  then  calculate 

m 

T  =  I  (e,  -  (29) 

1=1 

where  Q,  is  the  computed  discharge  at  the  /th  time. 

(5)  According  to  the  value  of  F,  the  values  of  K  and  SM  are  determined  for 
the  next  iteration  based  on  the  principle  of  the  simplex  method.  We  then 
return  to  step  (2).  The  details  of  computation  steps  are  shown  in  the  flow  chart 
(Fig.  1). 

Although  there  are  eight  parameters  in  the  rainfall-runoff  model,  only  two 
parameters,  representing  the  infiltation  process,  are  determined  by  iteration, 
and  for  this  purpose  it  is  convenient  to  use  the  simplex  method.  A  detailed 
discussion  of  this  method  has  been  given  by  Himmelblau  (1972). 

ESTIMATION  OF  PARAMETERS  OF  EXCESS  RAINFALL-RUNOFF  PROCESS 
FROM  MULTIPLE  EVENTS 


If  data  are  available  on  several  excess  rainfall-runoff  events,  the  average 
values  of  the  model  parameters  can  be  estimated  by  using  the  least-squares 
method  throughout.  We  assume  that  the  first  excess-rainfall  hyetograph  is 


^ii^:iAi  •  •  •  All  (30) 

and  the  corresponding  runoff  hydrograph  is 

QwQuQm  ■  ■  ■  Qn\  ■  ■  ■  Qm\  (31) 

The  second  excess-rainfall  hyetograph  is 

A; ^2; A: A’  •  •  •  Ai2  (32) 

and  corresponding  runoff  hydrograph  is 

Q12Q21QZ2Q42  ■  ■  •  Qm2  (33) 


According  to  the  least-squares  method,  the  average  model  parameters  are  as 
follows; 


^  =  [a^a2aybfjbih2]^ 


where 

rA,! 

B^B  =  [AjAj] 

1 

A2_ 

FT, I 

B^y  =  lAjAj] 

1 

y. 

-h  AJA2 

Aj  y,  -t-  Aj  y. 


(34) 

(35) 


(36) 
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read  data  Dt  (period) ,D=(rettention  capacity) 
I ( j ) =lnt ens ity ,  j  =number  of  period _ 


Initial  guess  parameters  K,  SM 
_ X<»=(K,SM)'^ _ _ 
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compute  new  estimates 
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V 

(40) 

APPLICATION 


The  proposed  model  was  applied  to  simulate  the  runoff  hydrograph  on  ten 
watersheds,  seven  of  which  are  located  in  the  loess  plateau  and  Sichan 
Province  in  China,  and  three  in  Texas,  USA,  respectively.  The  rainfall-runoff 
events  were  divided  into  two  groups.  One  group  was  utilized  to  calibrate 
the  model  and'  estimate  its  parameters,  and  the  other  was  used  for  model 
verification. 

Model  calibration 

The  model  chosen  for  simulating  the  excess  rainfall-runoff  process  has 
auto-regressive  order  three  and  moving  average  order  three.  Ten  rainfall- 
runoff  events  over  ten  watersheds  were  used  to  estimate  the  infiltration 
parameters.  In  the  estimation  of  parameters  by  the  simplex  method,  the  water 
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*  *  *  SIMUUTED  HYDROGRAPH 
+  -h  +  OBSERVED  HYDROGRAPH 

Fig.  2.  Comparison  of  observed  runoff  hydrograph  with  that  calculated  from  a  new  approach  for  the 
calibration  of  the  Dujia  event  (15-8-1966). 


balance  condition  had  to  be  satisfied.  Accordingly,  the  volume  of  the  runoff 
calculated  from  the  hydrograph  was  equal  to  that  of  the  excess  rainfall 
calculated  from  the  hyetograph.  The  values  of  the  parameters  obtained  for  the 
watersheds  are  shown  in  Table  1.  A  comparison  of  observed  and  calculated 
hydrographs  for  three  representative  sample  calibration  events,  ranging  from 
the  worst  calibration  to  the  best  calibration,  is  shown  in  Figs.  2-4.  For  further 
evaluation  of  the  calibration  performance  of  the  model,  the  goodness  of  fit  of 
a  computed  hydrograph  to  an  observed  hydrograph  was  determined  by  means 
of  the  following  four  criteria. 

(a)  Relative  squared  error  (RSE) 

m 

1  mo  -  Qiof 

RSE  =  - 

I  Q^-iO 

)=  I 

where  Q^{1)  is  the  computed  discharge  at  time  /.  Q(i)  is  the  observed  discharge 
at  time  i  and  m  is  the  number  of  discharge  ordinates  (RSE  represents  the 
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*  *  ♦  SIMUUMED  HYDROGRAPH 
+  +  +  OBSERVED  HYDROGRAPH 

Fig.  3.  Comparison  of  observed  runoff  hydrograph  with  that  calculated  from  a  new'  approach  for  the 
calibration  of  the  Xizhun  event  (28-8-1966). 


overall  shape  of  the  hydrograph,  and  if  RSE  =  0,  the  computed  hydrograph 
will  coincide  with  the  observed  hydrograph); 

(b)  Relative  error  in  estimated  peak  (£p) 


= 


\Q\  -  Q, 


(c)  Time  of  peak  discharge; 

(d)  Time  base  of  the  hydrograph. 

Table  2  gives  values  of  these  four  criteria  for  all  calibration  events.  The  table 
shows  that  the  RSE  generally  was  less  than  0.01 ;  for  only  one  event  it  exceeded 
0.02.  The  maximum  value  of  RSE  was  0.1166,  and  the  minimum  was  only 
0.00157.  The  difference  between  observed  and  computed  peak  discharge  was 
very  small.  The  average  value  of  £p  was  0.0239;  the  maximum  value  of  £p  was 
0.0645  and  the  minimum  was  only  0.003.  The  time  of  peak  discharge 
calculated  from  the  model  and  that  from  the  observed  data  were  the  same  for 
all  but  two  calibration  events.  The  computed  time  base  was  the  same  as  the 
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+  +  +  SIMULATED  HYDROGRAPH 
*  •  •  OBSERVED  HYDROGRAPH 

Fig.  4.  Comparison  of  observed  runoff  hydrograph  with  that  calculated  from  a  new  approach  for  the 
calibration  of  the  Lillie  event  (19-7-1979). 


TABLE  3 


Calibration  for  parameters  of  the  excess  rainfall-runoff  model  by  using  multiple  least-squares 
method 


Name  of 
basin 

Area  of 
basin  (km”) 

Parameters  of  excess-runoff  model 

^0 

A, 

f’: 

I 

Dujia 

96.1 

1.292 

-0.416 

-0.00771 

0.0353 

-0.0518 

0.1164 

0.968 

Xizhun 

49.0 

0.501 

0.I0I7 

-0.0978 

0.1024 

0.145 

0.140 

0.982 

Shejia 

4.26 

0.615 

-0.178 

0.0525 

0.183 

0.365 

-0.0537 

0.994 

Shanchuan 

21.0 

0.755 

-0.253 

0.161 

0.0245 

0.226 

0.0176 

1 .005 

Wuansun 

3.9 

1.074 

-0.546 

0.266 

0.108 

-0.207 

0.266 

0.962 

Fuxiq 

1.39 

1.318 

-0.495 

0.0552 

0.0373 

-0.00831 

0.0569 

0.964 

Shoal 

18.2 

0.222 

-  0.00623 

0.0999 

0.394 

0.250 

0.0942 

0.962 
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*  •  •  SIMULATED  HYDROGRAPH 
+  +  -E  OBSERVED  HYDROGRAPH 

Fig.  5.  Comparison  of  observed  runoff  hydrograph  with  that  calculated  from  a  new  approach  for  the 
verification  of  the  Xizhun  event  (15-8-1966). 


observed  value  for  each  calibration  event.  Therefore,  the  overall  features  of 
hydrograph  calculated  from  the  model  were  very  similar  to  those  observed. 

When  the  excess  rainfall  events  were  calculated,  two  or  more  than  two 
corresponding  runoff  events  in  the  same  basin  were  utilized  to  estimate  the 
average  parameters  of  the  excess  rainfall-runoff  model  by  using  the  multiple 
least-squares  method  as  discussed  above.  The  values  of  the  parameters  so 
obtained  are  shown  in  Table  3. 

Model  verification 

With  the  values  of  the  infiltration  parameters  obtained  by  calibration  for 
each  watershed,  the  excess  rainfall  hyetograph  was  computed  for  each  event. 
The  runoff  hydrograph  was  then  computed  for  each  excess-rainfall 
hyetograph.  using  the  runoff  model  parameters  obtained  for  multiple  events 
during  calibration.  To  test  the  model,  seven  rainfall-runoff  events  were  used. 
The  RSE  was  computed  for  the  computed  hydrograph  of  each  event.  The 
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•  •  *  SiMUL^iTED  HYDROGRAPH 
+  +  -i-  OBSERVED  HYDROGRAPH 

Fig.  6.  Comparison  of  observed  runoff  hydrograph  with  that  calculated  from  a  new  approach  for  the 
verihcalion  of  the  Shejia  event  (14-7-1966). 


results  of  computations  are  shown  in  Table  4,  and  for  three  representative 
sample  events  in  Figs.  5-7.  The  values  of  RSE  generally  were  less  than  0.1 1, 
with  the  average  being  0.0671.  The  maximum  value  of  RSE  was  0.1355,  and 
the  minimum  was  only  0.0147.  The  average  relative  error  of  the  estimated 
peaks  (£p)  was  0.145.  The  maximum  value  of  £p  was  0.336,  and  the  minimum 
was  0.0023,  and  the  peak  occurred  at  almost  the  correct  time  for  all  the 
verification  events.  The  time  base  of  hydrograph  calculated  from  the  model 
was  the  same  as  the  observed  value  for  each  verification  event. 

The  accuracy  of  the  model  depends  on  the  spatial  distribution  of  rainfall 
and  other  factors.  If  the  rainfall  was  distributed  uniformly  in  space,  the  model 
parameters  estimated  using  the  least-squares  method  will  approximate  the 
true  values.  For  one  event,  the  rainfall  center  may  be  located  in  the  lower 
reaches,  and  for  the  other  event  it  may  be  located  in  the  upper  reaches.  If  the 
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*  *  <  SIMULATED  HYDROGRAPH 
+  +  +  OBSERVED  HYDROGRAPH 

Fig.  7.  Comparison  of  observed  runoff  hydrograph  with  that  calculated  from  a  new  approach  for  the 
verification  of  the  Shanchua  event  (15-8-1966). 


parameters  estimated  from  the  one  event  are  used  to  predict  the  runoff 
hydrograph  of  the  other  event,  the  prediction  values  would  have  some  errors. 
Therefore,  multiple  excess  rainfall-runoff  events  were  applied  to  estimate  the 
average  values  of  the  model  parameters  which  were  then  used  to  predict  the 
runoff  hydrograph  of  other  events.  In  this  case,  the  results  would  be  better 
than  the  parameters  based  on  a  single  event.  Let  us  consider,  for  example,  the 
Xizhun  watershed.  The  runoff  hydrograph  of  15  August  1966.  was  calculated 
by  using  two  groups  of  parameters.  In  one  group,  the  parameters  were 
estimated  from  a  single  event  shown  in  Table  1  and  in  the  other  group  are  the 
average  parameters  Table  2.  RSE  =  0.1482  and  £p  =  0.31  for  a  single 
event,  and  RSE  =  0.0147  and  £p  =  0.0023  for  the  average  parameters. 
Therefore,  the  results  showed  that  the  parameters  based  on  multiple  events 
were  better  than  those  based  on  single  events. 
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SUMMARY  AND  CONCLUSIONS 

The  model  developed  in  this  study  combines  the  excess-rainfall  process  and 
the  runoff-concentration  process.  Parameters  of  these  processes  were  simul¬ 
taneously  estimated  by  minimizing  the  sum  of  squared  differences  between 
observed  and  computed  discharges,  subject  to  the  constraint  of  maintaining 
the  water  balance.  The  model  has  the  following  advantages; 

(1)  With  the  assumed  initial  values  of  K  and  SM,  the  excess  rainfall 
hyetograph  can  be  estimated  for  any  event  and  the  least-squares  method  can 
then  be  used  to  determine  the  parameters  of  the  excess  rainfall-runoff  process 
from  the  observed  hydrograph.  The  sum  of  squared  differences  between 
observed  and  computed  hydrographs  can  be  utilized  to  determine  the  values 
of  K  and  SM  for  the  next  iteration.  Thus,  parameter  estimation  is  simple  and 
easy. 

(2)  K  depends  on  the  soil  properties  and  can  be  taken  as  constant  for  a  fixed 
watershed.  SM  depends  on  the  capillary  potential  and  the  antecedent  soil 
moisture;  however,  SM  can  be  estimated,  based  on  laboratory  or  field 
observations. 

(3)  There  exists  an  enormous  volume  of  data  on  rainfall-runoff  events  for 
watersheds  in  various  countries.  These  data  can  be  analyzed  to  estimate  K  and 
SM,  and  to  establish  a  relationship  between  SM  and  antecedent  soil  moisture. 
The  prediction  of  excess  rainfall  can  be  carried  out  when  K  and  SM  are  given. 

(4)  Figures  2-4  and  5-7  show  that  observed  and  estimated  runoff 
hydrographs  are  in  good  agreement  for  calibration  as  well  as  verification 
events.  The  model  can  be  adapted  for  real-time  flow  forecasting. 
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ABSTRACT 

Wang,  G.-T..  Singh.  V.P.,  Guo.  C.  and  Huang,  K.X..  1991.  Discrete  linear  models  for  runoff  and  sediment 

discharge  from  the  Loess  Plateau  of  China.  J.  Hydrol..  127:  153-171. 

Discrete  linear  models  were  developed  for  estimating  runoff  and  sediment  discharge  hydrographs  from 
the  Loess  Plateau  of  China.  A  regression  equation  was  also  established  relating  runoff  rate  and  sediment 
discharge.  Tested  on  five  small  basins,  the  results  were  in  good  agreement  with  observations.  For  the 
discrete  linear  transfer  runoff  model,  the  values  of  the  integral  square  error  were  generally  less  than  1  “  o  for 
all  calibration  events,  and  around  10%  with  an  average  value  of  9.36%  for  all  verification  events.  For  the 
discrete  linear  transfer  sediment  model,  the  calibration  coefficient  of  determination  for  all  five  basins  w  as 
more  than  97%.  and  the  verification  R-  was  more  than  91%  with  an  average  of  94.3%. 


INTRODUCTION 

The  upper  and  middle  reaches  of  the  Yellow  River  in  China  flow  through 
a  wide  plateau  with  an  area  of  about  580  000  km*;  it  is  covered  with  loess  and 
red  loess.  The  Loess  Plateau  is  located  in  a  transition  zone  with  monsoons  in 
the  southeast  and  an  arid  climate  in  the  northwest.  It  borders  the  Qinghai- 
Tibet  Plateau  on  the  east  and  is  bounded  by  Mongolia-Xin  Plateau  on  the 
north.  In  the  winter,  the  region  is  controlled  by  high  pressure  systems  in 
Mongolia,  when  the  Arctic  cold  air  mass  moves  in,  causing  strong  breezes  and 
falling  temperatures.  In  the  summer,  subtropical  Pacific  air  masses  with  much 
more  moisture  reach  the  region  and  cause  rainfall.  Owing  to  differing  distance 
from  the  ocean,  the  east  to  west  variation  of  precipitation  is  large.  The 
distribution  of  precipitation  in  time  is  also  extremely  nonuniform.  For 
example,  the  precipitation  in  July.  August  and  September  is  about  70"  o  of  the 
annual  total.  The  precipitation  for  one  event  is  generally  about  10%  of  the 
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annual  total  but  as  much  as  43%  has  been  recorded.  The  maximum  precipi¬ 
tation  of  53.1mm  recorded  for  5  min  in  China  is  from  this  region.  The 
potential  evaporation  is  about  2500mm  year"'.  The  ground  water  is  more 
than  100  m  below  ground  level.  Therefore,  the  soil  stratum  seldom  develops 
saturated  conditions. 

Owing  to  the  particular  physiographic  characteristics  of  the  Loess  Plateau, 
the  distinguishing  features  of  runoff  from  rainfall  are  as  follows. 

(1)  Rivers,  with  drainage  areas  of  less  than  500  km‘,  are  nearly  dry  valleys. 
In  the  summer,  when  the  rainstorms  occur,  the  resulting  runoff  hydrographs 
are  characterized  by  short  durations  and  high  flood  peaks.  Generally,  there  is 
little  baseflow  to  be  separated  from  the  hydrographs.  Therefore,  discrete 
linear  transfer  function  models  can  be  suitable  for  simulating  the  hydrographs 
from  this  region. 

(2)  The  sediment  yield  for  the  basins  is  caused  by  storm  rainfall  and  runoff. 
Therefore,  the  relationship  between  sediment  and  runoff  discharge  is  very 
strong.  The  regression  equation  and  discrete  linear  transfer  models  are  proposed 
for  simulating  the  relation  of  sediment  discharge  to  the  runoff  hydrograph. 

A  large  amount  of  sediment  is  eroded  from  the  Loess  Plateau  and  enters  the 
Yellow  River.  Evaluation  of  this  yield  is  therefore  important  in  the  design  of 
soil  conservation  and  pollution  control  practices,  as  well  as  in  the  design  and 
management  of  dams,  canals  and  other  hydraulic  structures.  A  number  of 
studies  have  been  carried  out  to  investigate  soil  and  water  conservation 
aspects  of  the  Loess  Plateau  (Gong  and  Jiang,  1977).  Yin  and  Chen  (1989) 
analyzed  over  4000  observations  on  58  small  watersheds  with  basin  area  from 
0.193  to  329  km^  The  data,  covering  the  period  1954-1982,  included  21 
variables.  They  related  erosion  intensity  to  a  composite  index  of  basin  surface 
characteristics  (Yin  and  Chen,  1989). 

This  paper  develops  a  linear  discrete  model  to  simulate  runoff  and  sediment 
discharge  from  the  Loess  Plateau,  taking  into  account  its  physiographic 
characteristics. 

DISCRETE  LINEAR  TRANSFER  FUNCTION  MODELS  FOR  SMALL  BASINS  IN  THE 
LOESS  PLATEAU 

The  discrete,  linear  rainfall-runoff  model  given  by  Wang  and  Y u  ( 1986)  can 
be  expressed  as 

QU)  =  a,  (2(/  -/■)  +  •••  +  -  P)  +  +  bj{t  -  1) 

-l-  .  .  .  -P  b^I{t  —  (j)  ( 1 ) 

where  Q{t)  is  direct  runoff  discharge  at  time  /,  Q{i  -  1)  is  direct  runoff 
discharge  at  time  (/-/),/=  1.2 . p,  I{i)  is  excess  rainfall  at  time  t. 


RUNOFF  AND  SEDIMENT  DISCHARGE  FROM  THE  LOESS  PLATEAU 


155 


/(/  —  j),j  =  1,  2,  .  .  .  ,  (7,  is  excess  rainfall  at  time  (/  —  j),  and  Qj  and  bj  are 
coefficients. 

In  the  model  used  in  this  study,  the  excess  rainfall  (effective  rainfall)  is 
expressed  as  the  difference  between  rainfall  P{t)  and  losses  L{t)  over  the  same 
time  period,  that  is 

/(/)  =  P{t)  -  L(0  for  t  =  1,  2,  .  .  .  ,  5  (2) 

Inserting  eqn.  (2)  in  eqn.  (1)  yields 
Q(l)  =  bolPd)  -  L(l)] 

Q(2)  =  a,Q(l)  +  bolP(2)  -  L(2)]  +  b,lP(l)  -  W)] 


Q(s)  =  a,Q(s  -  1)  +  .  .  .  -h  a^Q(s  -  p)  +  bo[P(s)  -  L(s)] 
+  ...  +  b.lPis  -  q)-  Lis  -  q)] 


(3) 


Q(m)  =  a^Qim  -  1)  +  .  .  .  +  a^Qd  -  p)  J 

For  discrete  linear  models  represented  by  eqn.  (1),  the  parameters  can  easily 
be  determined  from  effective  rainfall  and  discharge  data  by  ordinary  least- 
squares.  correlation  analysis,  linear  programming  or  other  methods. 
However,  for  the  discrete  linear  models  given  by  eqn.  (3),  both  the  model 
parameters  a, .  o,. .  .  .  ,  and  the  time-variant  losses  L{t),  which 

are  not  known  beforehand,  must  be  estimated  simultaneously.  In  this  study, 
the  parameter  estimation  problem  is  formulated  as  a  nonlinear  minimization 
problem  as  follows. 

Let  the  residuals  between  Q'it)  and  estimated  Qit)  discharge  rates  be 
represented  by  e,  as 

e,  =  Q'ii)  -  0(0  for  /  =  1,2 . m  (4) 

substitution  of  eqn.  (3)  into  eqn.  (4)  yields  the  following  equation  for 
residuals: 


c,  =  b,[P{\)  -  KD]  -  Q(l) 

e.  =  a,Q{\)  +  h,[Pi2)  -  L(2)]  +  />,[/>(!)  -  1(1)]  -  Q(2) 


V  (5) 


-  1)  +  a.Qim  -  2)  -I-  ...  +  a^QUn  -  p)  -  Qim)  J 
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For  convenience,  let  the  unknown  a,  b  and  L{s)  values  in  eqn.  (5)  be  the 
components  of  a  column  vector  X 

X  =  [a,,  fl.,  .  .  .  ,  Qp,  Z>o,  ^1,  •  .  •  ,  L(l),  L(2),  .  .  .  ,  L(5)f  (6) 

where  T  denotes  the  transpose.  The  residuals  .  .  .  ,6>„ineqn.  (5)  are  then 
a  function  of  X.  The  optimal  estimate  of  A'  can  be  obtained  by  solving  the 
unconstrained  nonlinear  minimization  problem  with  the  following  objective 
function: 


F{X)  =  X  e]{X)  =  E\X)E{X)  (7) 

j  =  i 

where  E(X)  =  [e^{X),  e-<{X),  .  Taking  the  first  derivative  of 

eqn.  (7)  with  respect  to  X  yields 


where 


ce,(X) 

de,(X) 

de,(X) 

COi 

da^ 

dLis) 

ce.{X) 

ce.iX) 

de.(X) 

A(X)  = 

da^ 

cL{s) 

beJX) 

ceJX) 

rn, 

<5^2 

cLis) 

Taking  the  second  derivative  of  eqn.  (7)  with  respect  to  A’  yields 
dA'-  dA'-  dA^  dA'  ^  ) 


2d=£ 

d-A' 


Neglecting  the  second-order  term  gives 
d-£  ^  . 

dF  ^  ^  (9) 

Substituting  eqns.  (8)  and  (9)  into  the  Newton-Raphson  iterative  formula  yields 
dF(X) 

=  .V"  -  (A\4r'A^E  (101 

dA'- 


(A'  .4)-'A^E 
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When  eqn.  (10)  is  solved  iteratively,  two  problems  may  arise.  First,  the  matrix 
(A^ A)  is  sometimes  ill-conditioned  and  it  is  difficult  to  find  the  inverse. 
Second,  if  the  direction  of  search  and  the  gradient  of  the  vector  X  are 
approximately  orthogonal,  the  iterative  speed  is  very  slow.  In  order  to 
overcome  these  problems,  Marquardt  (1963)  suggested  a  revised  iterative 
method  which  changes  eqn.  (10)  to  the  following  form: 

;^(A+!)  ^  x^K)  _  [^T^  Wiy'A^E  (II) 

where  /  is  a  unit  matrix  and  W  is  constant.  When  W  =  0,  eqn.  (11)  reduces 
to  eqn.  (10).  With  W  sufficiently  large,  WI  can  overwhelm  A^A  and  the 
minimization  approaches  a  steepest  descent  search.  Therefore,  the  Marquardt 
method  can  be  seen  to  be  a  combination  of  the  Gauss-Newton  iterative 
method  and  steepest  descent.  This  method  forces  the  Hessian  matrix  to  be 
positive  definite  at  each  stage  of  minimization  and  ensures  that  the  estimate 
of  its  inverse  is  also  positive  definite.  The  algorithm  using  the  Marquardt 
method  is  described  by  Bard  (1970). 

RUNOFF  AND  SEDIMENT  DISCHARGE  FROM  SMALL  BASINS  OF  THE  LOESS 
PLATEAU 

The  sediment  yield  from  the  Loess  Plateau  basins  is  caused  by  storm 
rainfall  and  runoff.  The  greater  the  rainfall  intensity,  the  more  the  runoff;  and 
the  sediment  transporting  capacity  increases  with  runoff  discharge.  The 
influence  of  runoff  discharge  on  sediment  discharge  is  more  evident  than  that 
of  rainfall  and  its  intensity.  Two  methods  were  used  to  study  the  relationship 
between  runoff  discharge  and  sediment  discharge. 

Regression  method 

Gong  and  Jiang  (1977)  plotted  on  log-log  paper  the  relation  between 
observations  in  time  of  runoff  discharge  and  the  corresponding  sediment 
discharge  for  four  basins  of  the  Loess  Plateau.  They  concluded  that  for 
discharges  more  than  10m^s“',  the  relation  between  runoff  discharge  and 
sediment  discharge  was  very  close  to  a  straight  line,  but  was  erratic  for 
discharges  less  than  lOm'^s"'.  The  runoff  discharge  and  sediment  discharge 
data  from  six  small  basins  were  analyzed  in  this  study.  For  each  basin,  the 
regression  equation  was  derived  from  the  runoff  sediment  discharge  data  of 
more  than  two  events.  The  results  are  given  in  Table  1. 

Table  1  shows  that  the  exponent  of  each  regression  equation  approximates 
1.0.  with  the  average  exponent  index  for  the  six  basins  (of  area  from  0.18  to 
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187.0  km*)  being  0.089.  The  regression  coefficient  ranges  from  0.335  to  0.557. 
The  exponents  and  coefficients  do  not  depend  on  the  basin  area  because 
sediment  discharge  depends  mainly  on  storm  rainfall  and  runoff.  The  flow 
velocities  in  dry  valleys  are  high,  with  the  maximum  velocity  being  5.0- 
6.0  ms~'.  The  travel  time  for  the  basin  is  relatively  short  and  runoff  with 
sediment  reaches  the  outlet  quickly.  Field  investigations  showed  that  there 
was  little  sediment  deposited  in  the  dry  valley  during  the  flood  period. 

For  use  on  a  particular  basin,  the  runoff  discharge,  from  which  sediment 
discharge  rates  would  be  calculated  using  the  regression  equation,  can  be 
estimated  by  applying  the  discrete  linear  model  to  the  rainfall  data. 


Discrete  linear  model  for  sediment  discharge 

Caroni  et  al.  (1984)  applied  two  simple  stochastic  models  to  sediment  yield 
data  from  an  experimental  basin  in  the  U.S.A.  Wang  and  Yu  (1986)  proposed 
a  more  general  linear  discrete  transfer  function  model  of  the  type  introduced 
by  Box  and  Jenkins  (1976)  expressed  as 

Q{t)  =  a,Q{t  -  \)  ^  .  .  .  +  q^Qit  -  p)  +  b^Iit)  +  .  .  .  +b^I(t  -  q) 

(la) 

This  type  of  model  is  proposed  for  estimating  sediment  discharge  from  small 
basins  of  the  Loess  Plateau  in  China.  This  model  is  superior  to  regression 
equations  which  cannot  account  for  the  time-series  nature  of  rainfall,  runoff 
and  sediment  discharge  processes.  When  eqn.  (1)  is  applied  to  sediment 
discharge  rates  of  small  basins,  Q{t)  is  replaced  by  G(/),  and  /(/)  by  Q{t). 

G{t)  =  C,G(/  -  1)  +  .  .  .  +  C^Git  -  p)  +  W,Q{t)  +  ...+  W^Q(t  -  q) 

(12) 

where  C,,  . and  IFp,  ILj . are  parameters.  The  model 

parameters  can  be  estimated  from  observed  runoff  and  sediment  discharge 
data  using  ordinary  least-squares.  Let  e^(t)  be  the  difference  between  observed 
sediment  discharge  and  sediment  discharge  computed  from  eqn.  (12)  using  the 
observed  runoff  discharge  as  input 


eft)  =  G(t)  -  IC,G(/  -  1)  +  .  .  .  +  c^G(t  -  p)  +  ...  + 

+  ...  -h  U'^Qit  -  q)]  (13) 
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The  objective  function  for  estimating  the  parameters  C  and  fV  values  can  be 
expressed  as 


i  r G(o  -  i  QG{t  -  i)-  i  WjQit  -  j) 

/=i  L 


y=0 


(14) 


The  necessary  conditions  for  eqn.  (14)  to  have  an  optimal  solution  for  C  and 
W  values  are 


dF 

dC, 


o  I 

/=! 


C(<)  -  I  Cfl(l  -  r)  -  I  -  J) 


-12 


y=o 


dC, 


0 


for  /  =  U  2,  3,  .  .  .  ,  P 

m  m  P 

I  C(/)C((  -  k)  =  II  C,C((  -  i)G('  -  k) 

/=!  /=l  /  =  ! 

+  t  t  »',QU  -  j)G(l  -  k)  (15) 

;=l >=0 

letting 

m 

rccik  -  i)  =  I  G(t  -  i)G(t  -  k) 

/  =  I 
m 

^Qo(^  ^  J)  ^  X  j)G{t  —  k) 

i  =  I 

Equation  (15)  reduces  to 

''gg(^)  =  G^r^cik  —  1)  +  .  .  .  +  CprcciO)  + 

+  .  .  .  +  ~  ^)  (1^1 

In  the  same  manner,  let 


^  =  0  •  for;  =  0.1.2 . , 

which  yields 

—  G^rQQ{  —  l  +  1)  +  .  .  .  +  Cptf^Q^P  —  /)  +  W^^rQQ(l) 

+  .  .  .  +  W^rQQ{())  (17) 

The  model  parameter  can  be  estimated  from  a  subset  of  eqns.  (16)  and  (17) 
which  relate  the  model  parameters  to  the  correlation  functions  calculated 
from  observed  runoff  and  sediment  discharge  rate  series.  The  number  of 
eqns.  (16)  and  (17)  is  equal  to  the  number  of  model  parameters. 
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Setting  k  =  \,  1,  .  .  .  ,  P  and  /  =  0,  I,  2,  .  .  .  ,  ^  in  eqns.  (15)  and  (16), 
respectively,  yields  the  following  equations; 

'•cc(l)  =  C.raciO)  +  .  .  .  +  CprcdP  -  1)  +  W,rQ^{\) 

+  .  .  .  +  -  q) 

rccH)  =  C,rcc(l)  +  .  .  .  +  -  2)  +  WorQc(2} 

+  .  .  .  +  fV^rgc(2  —  q) 


(18) 


rociP)  =  C.r^ciP  -  1)  +  .  .  .  +  CprcdO)  +  W,rgc(P) 

+  .  .  .  +  ^^rgciP  -  q) 
and 

'■g£i(0)  =  C’i^ce(l)  +  •  •  •  +  CprQQ{P)  +  fV^rggiO) 

+  ...  +  J^^rgg(q) 

''Gq(-I)  =  C,rcg(0)  +  ...  +  Cpr^giP  -  1)  +  W^rgg^X) 

+  .  .  .  +  ^Qgiq  ~  1 ) 


rcQi-q)  =  C.rc^d  -  q)  +  ...  +  Cpr^giP  -  q)  +  W^rgg{q) 

+  .  .  .  +  fV^rgg(O) 

There  are  (p  +  ^  +  1)  parameters  in  the  same  number  of  linear  algebraic 
equations.  The  solution  of  eqn.  (18)  in  matrix  form  is 


B  =  A'D 


(19) 


A  = 


B  =  [C,,  C,  .  . 

.  ,  Cp,  Wo,  W,, 

. . . ,  w^Y 

^  Gc(0) 

■  •  raoiP  -  1) 

■  ^eG(i  ~  q) 

^GGd) 

•  •  r^ciP  -  2) 

^qg(2)  •  ■ 

■  '-qg(2  -  q) 

^GciP) 

•  •  ^gg(I^) 

^qg(P)  ■  • 

•  ^qg(p  ~  q) 

'"gq  ( 1 ) 

•••  rGg(P) 

^qq{^)  ■  • 

•  rgg(q) 

^G0(O) 

■  ■  •  ^GQiP  ~  1) 

rQQ{\)  .. 

■  fggiq  -  1) 

''Ge(i  -  q) 

•  •  •  ^gq(P  ~  q) 

'‘Qgiq)  ■  ■ 

-  D  =  [roG 

(1Kg(2)  .  .  . 

(PKq(O)  . . 

■  ''Ggi~q)V 
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APPLICATION  OF  MODELS 

The  two  discrete  linear  models  were  applied  to  rainfall-runoff  and  sediment 
discharge  data.  One  model  was  used  to  estimate  the  runoff  hydrograph  from 
observed  rainfall,  and  the  other  was  used  to  calculate  sediment  discharge  from 
computed  runoff. 

Rainf all-runoff-sediment  data 

The  data  used  in  this  study  were  collected  from  the  hydrologic  yearbook  for 
an  experimental  basin  designated  as  Chobagou  basin,  in  Zizhou  county, 
ShaanXi,  on  the  Loess  Plateau  of  China.  The  basin  area  is  187.0 km*.  There 
are  six  subbasins  with  areas  ranging  from  0.017  to  96.1  km*.  The  location  of 
these  subbasins  is  shown  in  Fig.  1.  Rainfall  at  28  gauges  within  the  basin,  and 
runoff  and  sediment  concentration  at  their  outlets  are  recorded  in  data  files  at 
unequal  time  intervals.  A  5  min  interval  was  chosen  for  all  the  subbasins 
except  Twanshan  and  Twanshan  (No.  9)  (1  min  interval)  in  order  to  get  a 
sufficiently  detailed  resolution  of  data.  The  location  of  gauging  stations  for 
rainfall  and  runoff  discharge  is  also  shown  in  Fig.  1.  The  mean  areal  rainfall 
was  estimated  as  the  arithmetic  mean  of  the  station  records  in  the  neighbor¬ 
hood  of  the  centroid  of  the  subbasin. 
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Thirteen  rainfall-runoff-sediment  events  were  selected  to  assess  the 
adequacy  of  the  discrete  linear  models  for  simulating  runoff  and  sediment 
hydrographs.  The  data  were  divided  into  two  sets;  a  calibration  set,  and  a 
prediction  set.  The  data  in  the  calibration  set  consisted  of  eight  events  which 
were  used  for  parameter  estimation.  The  data  in  the  prediction  set  consisted 
of  five  events  which  were  used  for  model  verification. 

Calibration  of  runoff  hydrograph  model 

The  model  chosen  has  autoregressive  order  three  and  moving  average  order 
three.  Eight  rainfall-runoff  events  over  five  basins  were  used  to  estimate  the 
model  parameters  and  losses.  There  were  thus  five  model  parameters  and 
losses  L(l),  L(2),  .  .  .  ,  L{s),  to  be  determined  in  the  model.  The  number  of 
interval  losses  depends  on  that  of  rainfall  intervals.  For  minimizing  the  residual 
sum  of  squares  of  differences  between  observed  and  estimated  discharges, 
Marquardt’s  (1963)  iterative  algorithm  for  computing  nonlinear  least-squares 
solutions  was  used.  The  algorithm  needs  some  starting  values  for  the  parameters 
and  criteria  for  convergence.  The  goodness  of  fit  of  the  computed  hydrograph  to 
the  observed  hydrograph  was  estimated  by  means  of  the  following  three  criteria; 

(1)  Integral  square  error  (ISE) 

IQV)  -  G(()]= 

ISE  =  -  X  100% 

I  2(0 

/=  I 

(2)  Relative  error  in  estimated  peak  (EP) 

EP  =  ~  X  100% 

(3)  Time  of  peak  discharge. 

Table  2  gives  the  three  criteria  for  all  calibration  runs. 

Table  2  shows  good  agreement  between  observed  and  simulated  hydro¬ 
graphs  for  all  calibration  events.  The  values  of  ISE  generally  are  less  than 
1.0%;  for  only  two  events  did  ISE  exceed  5.0%.  The  maximum  value  of  ISE 
is  6.3%.  and  the  minimum  is  only  0.03%.  The  differences  between  observed 
and  computed  peak  discharges  are  very  small  and  the  times  of  peak  discharge 
were  also  well  reproduced  for  all  calibration  events.  The  model  parameters  for 
the  five  small  basins  are  given  in  Table  3. 

Verification  of  runoff  hydrograph  model 

The  values  of  the  parameters  from  calibration  events  for  each  basin  were 
applied  to  calculate  runoff  hydrographs  using  the  excess  rainfall,  which  was 
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TABLE  2 


Comparison  of  computed  and  observed  hydrographs  for  calibration  events  of  the  Loess  Plateau 
region 


Basin 

Area 

(km-) 

Event 

Peak  discharge 
(m^s"') 

EP 

(%) 

Peak  time 

ISE 

Observed 

.  r/o) 

Computed 

Observed 

Computed 

Dujiago 

96.1 

668028 

419 

423 

0.9 

19:25 

19:25 

5.63 

660815 

606 

607.3 

0.3 

19:05 

19:05 

0.65 

Xizhuang 

49 

660717 

280 

280 

0.0 

18:35 

18:35 

1.90 

660627 

16.0 

15.9 

0.4 

15:05 

15:05 

6.30 

Shanchuan 

21 

660717 

40.8 

40.8 

0.0 

18:45 

18:45 

0.63 

660828 

109.7 

109.8 

0.09 

19:20 

19:20 

0.80 

Shejia 

4.26 

660815 

44.1 

44.1 

0.0 

19:25 

19:25 

0.03 

Twanshan 

0.18 

660717 

2.08 

2.07 

0.3 

18:39 

18:39 

0.48 

estimated  from  the  rainfall  and  loss  processes  obtained  from  calibration 
events.  In  order  to  assess  the  accuracy  of  the  discrete  linear  model,  five 
rainfall-runoff  events  were  used.  The  results  are  given  in  Table  4. 

Table  4  shows  the  results  of  applying  the  discrete  linear  hydrograph 
model  for  all  verification  events.  The  values  of  ISE  generally  are  around 
10%  with  the  average  being  9.36%.  The  maximum  value  of  ISE  is  12.9%.  and 
the  minimum  is  only  4.4%.  The  average  relative  error  in  estimated  peak  (EP) 
is  18.92%.  The  maximum  value  of  EP  is  33.4%,  and  the  minimum  is  only 
3.3%.  and  the  peak  occurred  at  almost  the  correct  time  for  all  the  verification 
events. 


TABLE  3 


Calibration  results  for  runoff  hydrograph  model 


Basin 

Area 

(km-) 

Parameters 

^2 

^0 

b: 

Dujiago 

96.1 

1.121 

-0.252 

-0.510 

Xizhuang 

49.0 

0.3674 

0.0962 

0.0294 

0.1608 

Shanchuan 

21.0 

1.0387 

-0.4346 

0.1574 

0.0062 

-0.0341 

0.2636 

Shejia 

4.26 

0.6127 

-0.3316  , 

0.0637 

0.0004 

0.4822 

0.1696 

Twanshan 

0.18 

1.341 

-0.642 

0.076 

0.282 

-0,091 

0.023 
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TABLE  4 


Verification  results  for  runoff  hydrograph  model 


Basin 

Area 

(km-) 

Event 

Peak  discharge 
(m^s"') 

EP 

(%) 

Peak  time 

ISE 

(%) 

Observed 

Computed 

Observed 

Computed 

Dujiago 

96.1 

660815 

621.2 

641.7 

3.3 

19:30 

19:30 

7.0 

Xizhuang 

49.0 

660828 

446.6 

297.4 

33.4 

19:00 

19:00 

10.7 

Shanchuan 

21.0 

660815 

354.8 

274.7 

22.6 

19:20 

19:15 

12.9 

Shejia 

4.26 

660717 

33.1 

38.7 

16.9 

18:45 

18:40 

11.8 

Twanshan 

0.18 

660815 

6.97 

5.69 

18.4 

19:17 

19:16 

4.4 

Calibration  of  sediment  discharge  models 


The  equations  for  the  regression  of  sediment  discharge  on  runoff  discharge 
are  given  in  Table  1 ,  which  can  be  used  for  estimating  sediment  discharge  from 
computed  runoff  discharge. 

For  application  of  the  discrete  linear  model,  sediment  discharge  data  from 
five  small  basins  were  utilized  for  estimation  of  parameters.  The  model  chosen 
has  autoregression  and  moving  averages  each  of  order  two.  That  is 

G{t)  =  C,G(t  -  1)  +  C,G{t  -  2)  +  lV,Qit)  +  fV,Q(t  -  1)  (20) 

The  parameters  of  eqn.  (20)  were  estimated  from  one  event  for  each  basin  by 
least-squares.  The  residual  square  error  (RSE)  was  used  to  compute  the 
proportion  of  variance  explained  by  the  coefficient  of  determination 

(RSE/5,)= 

where  S,  is  the  output  standard  error  of  estimate.  Five  events  were  selected  for 
the  calibration  for  five  tested  basins.  The  model  parameters  and  calibration 
coefficient  of  determination  for  the  five  basins  are  given  in  Table  5.  as  are 
the  results  of  calibration.  The  good  agreement  expressed  by  high  values  of  the 
coefficient  of  determination  is  compatible  with  a  low  range  of  residuals. 


Prediction  of  sediment  discharge 

The  values  of  the  parameters  from  the  calibration  events  for  each  basin 
were  applied  to  calculate  sediment  discharge  using  the  estimated  runoff 
discharge.  For  the  regression  equation,  the  coefficient  and  exponent  were 
determined  from  the  data  of  more  than  two  events.  In  order  to  assess  the 
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TABLE  5 


The  calibration  results  for  sediment  discharge  rates 


Basin 

Area 

(km-) 

Event 

Parameter 

PC' 

C, 

C, 

% 

W, 

Dujiago 

96.1 

660828 

0.433 

-0.0703 

0.275 

0.01224 

0.9971 

Xizhuang 

49.0 

660717 

-0.175 

0.1596 

0.2893 

0.0632 

0.9715 

Shanchuan 

21.0 

660828 

-0.676 

0.0717 

0.2156 

0.4250 

0.9975 

Shejia 

4.26 

660815 

0.4908 

-0.079 

0.3287 

-0.1421 

0.9925 

Twanshan 

0.18 

660717 

0.5705 

-0.1333 

0.3805 

-0.1518 

0.981 

accuracy  of  the  two  models,  five  rainfall-runoff-sediment  discharge  events 
were  used.  The  results  are  given  in  Table  6. 

Table  6  shows  that  for  verification  events,  both  models  gave  good  results. 
For  the  discrete  linear  function  model,  the  value  of  R~  is  more  than  91.0% 
with  an  average  value  of  94.26%.  The  maximum  is  96.4%,  and  the  minimum 
is  91.0%.  For  regression  equations,  the  agreement  expressed  by  high  values  of 
the  coefficient  of  determination"  is  comparable  with  that  of  discrete  linear 
transfer  model. 

TABLE  6 

Verification  results  for  sediment  discharge  rates 


Basin  Area  Event  Peak  sediment  discharge  EP  R‘ 

(km-)  (m^s"') 

-  DLFM  RE  DLFM  RE 

Observed  Computed 


DLFM  RE 


Dujiago 

96.1 

660815 

280.0 

262.7 

261.5 

6.2 

6.6 

96.4 

95.4 

Xizhuang 

49.0 

660828 

145.0 

96.5 

122.2 

33.4 

15.7 

96.4 

98.2 

Shanchuan 

21.0 

660815 

136.8 

108.3 

97.7 

20.8 

28.6 

91.1 

93.4 

Shejia 

4.26 

660717 

12.6 

13.1 

13.0 

4.0 

3.2 

91.0 

91.0 

Twanshan 

0.18 

660815 

3.1 

2.35 

1.97 

24.2 

36.5 

96.4 

93.6 

DISCUSSION  OF  RESULTS 

Discrete  linear  transfer  function  models  for  runoff  discharge  were  used  to 

estimate  the  parameters  a,,  fl; . Op,  .  .  .  ,  and  L(l),  L(2) . 

L{s).  and  these  parameters  were  used  to  calibrate  the  model.  The  calibration 
results  are  in  excellent  agreement  with  observations.  Strictly  speaking,  these 
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parameters  cannot  directly  be  used  for  another  rainfall-runoff  event  for  the 
following  reasons. 

(1)  The  number  of  rainfall  periods  in  the  events  used  for  calibration  is  not 
equal  to  those  of  the  events  used  for  verification. 

(2)  The  infiltration  intensity  is  related  to  the  rainfall  intensity  on  the  Loess 
Plateau  based  on  the  analysis  of  the  data  from  observed  storm  runoff  in  small 
basins  and  artificial  rainfall  tests  in  the  field  (Liu  and  Wang,  1980).  The  loss 
of  rainfall  is  mainly  related  to  rainfall  intensity  {a)  and  soil  flow  properties. 
The  average  loss  rate  ^  can  be  plotted  against  the  average  intensity  {a)  as 

=  R’  (f 

where  R  is  the  coefficient  of  loss  related  to  the  kind  of  soil  or  soil  flow 
properties,  and  r  is  the  index  of  loss. 

(3)  The  condition  of  rainfall-runoff  events  for  calibration  is  not  the  same 
as  that  of  rainfall-runoff  events  for  verification. 

(4)  The  soil  moisture  of  rainfall-runoff  events  for  calibration  is  not  the 
same  as  that  of  rainfall-runoff  events  for  verification. 

Therefore,  this  method  in  practice  cannot  directly  be  used  for  prediction  of 
runoff  discharge.  However,  the  method  can  be  applied  to  analyze  the  loss 
process  from  the  hydrograph  for  calibration  events,  which  can  then  be  used 
to  estimate  the  loss  process  for  verification  events  based  on  the  rainfall 
intensity  and  soil  moisture,  and  excess  rainfall  can  be  calculated  from  the 
rainfall  and  loss  process.  In  this  way,  the  model  parameters  from  the  cali¬ 
bration  event  can  be  utilized  to  calculate  the  runoff  hydrograph  using  the 
excess  rainfall  for  the  verification  event. 

The  rainfall-runoff-sediment  relationship  is  time-dependent.  Consequently, 
the  model  parameters  for  a  particular  event  will  differ  from  the  average 
parameter  values.  In  order  to  overcome  this  drawback,  simultaneous  esti¬ 
mation  of  model  parameters  for  more  than  two  isolated  events  was  performed 
to  achieve  better  parameter  values. 

With  a  sufficient  number  of  calibration  events,  the  dependency  of  the 
models  as  well  as  parameter  values  on  hydrological,  meteorological,  veg- 
etational  or  other  characteristics  of  the  basin  can  be  investigated. 

SUMMARY  AND  CONCLUSIONS 

A  discrete  linear  model  was  employed  for  simulating  runoff  and  sediment 
hydrograph  for  five  small  basins  in  the  Loess  Plateau  of  China.  The  model 
parameters  were  estimated  using  a  nonlinear  least-squares  method  in  con¬ 
junction  with  Marquardt’s  (1963)  iterative  algorithm  from  the  calibration 
events.  The  results  showed  good  agreement  between  observed  and  estimated 
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Fig.  2.  Comparison  of  observed  runoff  and  sediment  discharge  hydrographs  with  those  calculated  from 
discrete  linear  models  for  a  calibration  event  on  17  July,  1966,  in  Shanchuan  basin. 


Observed  runoff 


Cdlcuiaced  runoff 
Observed  sediment 
Calculated  runoff 


Fig.  3.  Comparison  of  observed  runoff  and  sediment  discharge  hydrographs  with  those  calculated  from 
discrete  linear  models  for  the  verification  event  of  15  August,  1966  in  Dujiago  basin. 
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Fig.  4.  Comparison  of  observed  runoff  and  sediment  discharge  hydrographs  with  those  calculated  from 
discrete  linear  models  for  the  verification  event  of  15  August,  1966  in  Twanshan  basin. 

hydrographs  for  the  five  tested  basins.  The  values  of  ISE  generally  were  less 
than  1.0%.  The  coefficient  and  exponent  of  the  regression  equation  were 
determined  using  the  data  from  more  than  two  events,  and  the  parameters  of 
the  discrete  linear  model  for  sediment  were  estimated  using  ordinary  least- 
squares  in  the  calibration  events.  The  calibration  coefficient  of  determination 
R‘  for  five  small  basins  was  more  than  97%.  Figure  2  shows  the  relationship 
between  observed  and  estimated  runoff  and  sediment  discharge  for  a  sample 
event.  The  agreement  between  computed  and  observed  values  lends  credence 
to  the  models. 

The  results  from  the  verification  events  were  not  as  good  as  those  from  the 
calibration  events,  but  the  prediction  accuracy  of  the  discrete  linear  models 
both  for  runoff  and  sediment  discharge  was  satisfactory. 

For  the  discrete  linear  transfer  model  of  the  runoff  hydrograph,  the  values 
of  ISE  were  around  10%  with  the  average  value  being  9.36%.  For  the 
sediment  models,  both  discrete  linear  model  and  regression  equation  gave 
good  results:  the  value  of  R‘  for  both  were  more  than  91.0%.  The  regression 
equation,  with  the  average  value  of  R'  being  94.32%.  was  a  little  better  than 
the  discrete  linear  model,  where  the  average  value  of  R~  was  94.26%.  However. 
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for  relative  error  of  peak  sediment  discharge,  the  discrete  linear  model  gave 
better  results;  the  average  value  of  EP  was  17.72%  compared  with  18.12%  for 
the  regression  equation. 

Figures  3  and  4  show  the  relationship  between  observed  and  estimated 
runoff  and  sediment  for  sample  events  selected  from  verification.  The  good 
agreement  between  observed  and  predicted  runoff  and  sediment  discharge  so 
obtained  support  the  adaptability  of  these  models  to  real-time  forecasting 
schemes  in  the  Loess  Plateau  in  China. 
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APPENDIX:  NOTATION 

parameters  for  model  (1) 


parameters  for  model  (12) 

• 

e,„  the  difference  between  the  observed  and  estimated 
runoff  discharges 


ho. 

C,.  G. 

m;.  h; 

e,.  . 
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^g(0 

the  difference  between  the  observed  and  estimated 
sediment  discharge  rates 

Git) 

observed  sediment  discharge  rate 

GV) 

estimated  sediment  discharge  rate 

no 

excess  rainfall  intensity  (m^s“') 

Us) 

the  losses  of  rainfall 

m 

the  number  of  runoff  discharge  periods 

P,Cj 

the  order  of  models  (1),  (3)  and  (12) 

Pit) 

rainfall  intensity  (m^s“') 

Qp 

observed  peak  discharge 

Q'p 

modeled  peak  discharge 

Qit) 

observed  runoff  discharge  (m^s'') 

Q'it) 

estimated  runoff  discharge  (m^s“') 

S 

the  number  of  non-zero  rainfall  periods 

S.. 

estimated  standard  error  of  the  sediment  discharge 

rates 

Derivation  of  Infiltration  Equation  Using 
Systems  Approach 

By  V.  P.  Singhs'  Member^  ASCE,  and  F.  X.  Yu" 


Abstract:  A  general  infiltration  model  is  derived  using  a  systems  approach.  The 
models  of  Honon,  Kostiakov,  Overton.  Green  and  Ampt.  and  Philip  are  some  of 
the  example  models  that  are  shown  as  special  cases  of  the  general  model.  An 
equivalence  between  the  Green- Ampt  model  and  the  Philip  two-term  model  is  shown. 
The  general  model  also  provides  a  solution  for  the  Holian  model  expressing  in¬ 
filtration  as  a  function  of  time.  This  solution  of  the  Holtan  model  has  not  been 
reported  in  the  literature.  A  first-order  analysis  is  performed  to  quantify  the  un¬ 
certainty  involved  with  the  generalized  model.  The  general  infiltration  model  con¬ 
tains  five  parameters.  Two  of  the  parameters  are  physically  based  and  can  therefore 
be  estimated  from  the  knowledge  of  soil  properties,  antecedent  soil  moisture  con¬ 
ditions.  and  infiltration  measurements:  the  remaining  three  can  be  determined  using 
the  least  squares  method.  The  model  is  verified  using  ten  observed  infiltration  data 
sets.  Agreement  between  observed  and  computed  infiltration  is  quite  good. 


Introduction 

A  multitude  of  infiltration  models  used  in  applied  hydrology  and  soil  sci¬ 
ence  exists.  Some  of  these  models  are  theoretic^ly  based  (Dooge  1973:  Philip 
1957,  1969:  Green  and  Ampt  1911),  while  others  are  empirical  (Kostiakov 
1932:  Honon  1938:  Holtan  1961;  Overton  1964).  Some  of  the  empirical 
models  are  quite  popular  and  frequently  used  in  various  water  resources  ap¬ 
plications.  The  reason  for  their  popularity  is  that  they  are  simple  and  yield 
satisfactory  results  in  some  cases.  Since  the  empirical  models  are  more  or 
less  based  on  experimental  observations,  they  represent  the  overall  infiltra¬ 
tion  process. 

The  infiltration  process  can  be  represented  by  a  systems  approach,  which 
uses  an  absorber  or  a  network  of  absorbers.  Infiltration  constitutes  input  to 
the  absorber.  An  absorber  can  be  defined  by  a  relation  between  its  moisture 
content  and  the  rate  of  infiltration  to  it.  In  this  study,  a  general  definition 
of  the  absorber  is  suggested,  and  a  general  infiltration  model  is  denved  by 
coupling  this  general  relation  with  the  spatially  lumped  form  of  continuity 
equation  to  be  satisfied  by  the  absorber.  The  various  well-known  infiltration 
models,  such  as  the  models  of  Horton  (1938),  Kostiakov  (1932T,  Overton 
(1964),  Holtan  (1961),  Green  and  Ampt  (1911),  and  Philip  (1957)  can  be 
shown  as  special  cases  of  this  general  model.  Therefore,  the  systems  ap¬ 
proach  presents  a  unified  framework  for  all  of  these  models  and  shows  that 
these  models  result  from  different  definitions  of  the  absorber.  As  a  conse¬ 
quence.  connections  between  some  of  these  models  as  well  as  their  param¬ 
eters  can  be  established.  This  may  help  estimate  model  parameters.  In  par- 
ticular.  an  equivalence  between  the  Green-Ampt  model  and  the  Philip  two- 
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term  model  can  be  derived.  This  approach  also  leads  to  an  explicit  solution 
of  the  Holtan  model,  which  does  not  appear  to  have  been  reported  in  the 
literature. 

The  objective  of  this  study  is  to  use  the  systems  approach  to  derive  a 
generalized  infiltration  model,  show  the  various  infiltration  equations  re¬ 
ported  in  the  literature  as  special  cases  of  the  generalized  model,  show  con¬ 
nections  among  some  of  these  infiltration  models,  and  derive  an  explicit 
solution  of  the  Holtan  model.  The  study  also  presents  a  first-order  analysis 
of  uncertainty  and  develops  a  simple  procedure  for  estimating  model  param¬ 
eters.  Field  data  are  finally  used  to  verify  the  proposed  infiltration  model. 

Governing  Equations 

We  consider  a  column  of  soil  matrix  with  unit  area  of  vertical  infiltration 
as  shown  in  Fig.  1.  The  initial  storage  space  available  in  the  soil  column  is 
denoted  as  Sq,  If  the  soil  is  initially  dry.  then  5n  equals  effective  porosity 
multiplied  by  the  total  volume  of  the  column.  At  time  r  =  0,  a  thin  layer 
of  water  starts  to  cover  the  upper  surface  of  the  column.  At  any  time  t  the 
infiltration  rate  is  denoted  as  /(r),  seepage  rate  as  /,(r),  and  the  potential 
water  storage  space  as  S(t).  The  spatially  lumped  continuity  equation  for  the 
soil  column  can  be  written  in  integral  form  as 

5(/)  =  5o  -  I  findt  +  I  f,(r)dT . (D 

Jo  Jo 

The  differential  form  of  continuity  equation  can  be  obtained  by  differen¬ 
tiating  Eq.  1  with  respect  to  r  as 

dS(r) 

—  =fsU)  -  fir) . (2) 


FIG.  1.  Systems  Representation  of  Soil  Column  by  Absorber  with  Varying  Ver¬ 
tical  Infiltration  and  Seepage  Rate 
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The  amount  of  water  stored  in  the  soil  column  from  r  =  0  to  any  time 
1  is 


mt)  =  So-su)= 


If/jU)  =  0  in  Eq.  2  (Dooge  1973),  then 


dW 

—  =m 

dt 


(3) 


(4) 


Eq.  2  contains  three  unknowns:  S,f,  and  f^.  To  solve  this  equation,  two 
relationships  among  the  variables,  5(r),/(r),  and  fs(t),  are  required.  By  anal¬ 
ogy  to  rainfall-runoff  modeling  (Kulandaiswamy  1964),  a  possible  general 
relationship  among  these  variables  may  be  expressed  as 


dt  dt 


dJs 


A/  >  0. 


(5) 


where  a„  and  —  coefficients,  which  may  be  functions  of  time  /.  infiltration 
rate /(O,  and  seepage  rate  fs{t)\  and  M  =  some  integer. 

It  is  not  clear  what  the  form  of  the  third  relationship  should  be.  Of  course, 
one  can  arbitrarily  choose  some  relationship,  linear  or  nonlinear,  to  represent 
the  relation  among  the  three  variables  and  then  determine  their  coefficients 
by  experimental  data.  However,  the  resulting  model  (Eqs,  2.  5,  and  a  third 
relationship)  may  be  too  general  and  too  complicated  to  be  of  any  practical 
value.  Thus,  it  may  be  preferable  to  examine  the  existing  infiltration  models 
and  then  find  an  underlying  relationship  among  the  potential  storage,  infil¬ 
tration,  and  seepage  variables. 


Analysis  of  Some  Infiltration  Models 


Overton  (1964)  has  shown  that  several  infiltration  models  are  based  on  a 
relationship  between  infiltration  rate  (or  excess  infiltration  rate)  and  the  vol¬ 
ume  of  either  actual  or  potential  infiltration  (or  excess  infiltration).  For  ex¬ 
ample,  a  fonrb  of  the  Green  and  Ampt  model  (191 1)  can  be  derived  bv  as¬ 
suming  that 


/(O  =  —  +/. 
Fit) 


dFit)  a 

or  - = -  -h  L 

dt  Fit) 


(6) 


where  F  volume  of  infiltration  —  /o  fit)dt:  /.  =  the  ultimate  infiltration 
rate;  and  a  =  constant  of  proportionality.  Integration  of  Eq,  6  yields  a  Green 
and  Ampt  type  model; 


t 


(7) 


Similarly,  the  Kosliakov.  Philip  two-term.  Horton,  and  Overton  models  can 
be  derived.  Table  I  summarizes  these  models  and  the  postulates  on  which 
they  are  based.  In  the  approach  followed  by  Overton  (1964).  no  consider¬ 
ation  is  given  to  the  continuity  Eqs.  I  and  2. 
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TABLE  1.  Summary  of  Some  Infiltration  Models  Derived  by  Postulating  Reia' 
tionship  Between  Rate  of  Infiltration  and  Volume  of  Infiltration 


Assumed  relation 
(1) 

Ultimate 
(constant) 
rate  of 
infiltration, 

/ 

(2) 

Resulting  equation 
(3) 

Analogous 

to 

(4) 

fit)  =  a/\F{t)\ 

0 

F  -  \  Zai  or  y'  =  \  2/ a 

Kostiakov 

a  =  constant  of  proportional ity 

model 

Jin 

f. 

/  =  (l//)lF  -  (£///)ln(l  -  lF/ui/f.)]\) 

Green-Ampi 

=  a/lFil)\ 

a  —  constant  or  proponionality 

model 

/(')  - 

/ 

F  =  ^  \  2a!  or /  =  f]  +  Vu75  /*''* 

Philip  two- 

=  a/[F(l)  -//] 

term 

model 

/U'  -  / 

f 

f  -  f  ^  ifu  ~  t\  )  exp  <  -at) 

Honon 

= -  f. 

F  =  ft  +  [{ft  -  J[)/a\[\  -  exp(-ii/l] 

model 

-  uino  -  /  /) 

a  =  parameter,  ft  -  initial  infiltration  rate 

/«»  -  f. 

L 

/  =  ft  sec*((/,  -  f)  \  a/\\ 

Ovenon 

=  u[J,l  - 

=  ( 1  /\  at )  tan  '[F,  (a/f,  )*'"] 

model 

F^  =  ultimate  volume  of  infiltration 

These  models  plus  some  others  can  be  derived  using  the  systems  approach 
of  Eqs.  2  and  5.  According  to  this  approach,  the  potential  storage  S(f)  and 
the  seepage  rate  /,(r)  are  postulated  for  each  model.  The  relationships  ex¬ 
pressing  these  postulates  are  then  substituted  into  Eq.  2  to  obtain  the  cor¬ 
responding  infiltration  model. 

A  summary  of  eight  popular  infiltration  models  is  given  in  Table  2.  These 
models  are  presented  in  two  ways:  ( 1 )  The  usually  reported  form  in  literature 
(column  6);  and  (2)  the  form  derived  by  systems  approach  (column  5).  The 
potential  storage  function  S(f)  and  seepage  rate  function  f,(f)  for  each  model 
arc  also  given  in  Table  2:  a.  b,  A.  Sp.  Wo.  and  n  are  constants;  /  and  K 
denote  steady-state  infiltration  rate;  h  is  the  depth  of  ponded  water  on  the 
soil  surface  (a  constant);  6,  is  volumetric  water  content  in  transmission  zone; 
and  S,  is  an  effective  matrix  suction  head  at  the  wetting  front. 

It  may  be  pertinent  to  comment  briefly  on  the  Zhao  model,  which  is  com¬ 
monly  used  in  China.  2^ao  (1981)  hypothesized  the  infiltration  system  to 
be  comprised  of  an  absorber  and  a  regulator.  The  absorber  was  assumed  to 
be  inversely  proportional  to  the  amount  of  water  stored  in  a  finite  soil  col¬ 
umn.  and  the  regulator  to  be  directly  proportional  to  the  same  amount  of 
water.  This  hypothesis  leads  to  the  Zhao  model  as  shown  in  Table  2. 

From  column  2  of  Table  2.  it  is  seen  that,  except  for  Zhao's  model  ( 1981 ). 
/,())  is  either  equal  to  zero  (the  Kostiakov  model)  or  equal  to  the  steady- 
state  infiltration  rate.  If  we  compare  these  models  with  a  typical  infiltration 
curve  shown  in  Fig.  2.  then  /(r)  actually  represents  the  steady-state  How 
rate,  which  is  either  equal  to  a  constant  or  zero.  From  this  point  of  view. 
Zhao's  model  may  not  be  appropriate  because  the  steady-state  infiltration 
rate  may  not.  in  general,  be  proportional  to  the  amount  of  water  stored  in 
the  soil  column.  Therefore,  a  generalized  relationship  sought  between  fit) 
and  Sit)  should  not  conform  to  this  model.  .As  seen  from  Fig.  2.  the  quantity 
[/(/)  - /,(t)l  actually  represents  the  nonsteady  infiltration  rate,  and  this  func¬ 
tion.  according  to  column  3  of  Table  2.  is  either  directly  proponional  to  the 


840 


potential  storage  5(f)  (Horton  1938:  Overton  1964;  Holtan  1961)  or  inversely 
proportional  to  the  amount  of  water  stored  in  the  soil  column.  WU)  =  [5„ 
-  S(f)]  (Philip  1957;  Kostiakov  1932;  Green-Ampt  1911). 


Generalized  iNnuTRATiON  Model 


Based  on  this  analysis,  we  might  propose  a  generalized  infiltration  model 
of  the  following  form  of  which  all  seven  of  these  popular  models  are  special 
cases; 


/(/)  -m  = 


a|5(f)r 

[So  -  SU)]" 


(8) 


where  a.  m.  and  n  =  positive  real  constants.  According  to  Eq.  8. /(f)  — »  ^ 
at  f  =  0,  and/(r)  —*f,U)  at  /  =  ==.  The  general  form  of fjt)  may  be  proposed 

as 


where/  =  a  constant.  Eqs.  2.  8.  and  9  constitute  the  generalized  infiltration 
model. 

Substituting  Eq.  8  into  Eq.  2  yields 


dS(n  -al5(f)]'" 

dt  ~  15,  -  Sif)]" . 

Eq.  8  can  be  expressed  with  the  use  of  Eq.  9  as 
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(11) 


fit)  =f.  + 


a[5(/)]'" 
[So  -  S(t)r 


Integrating  Eq.  10 
[5o  -  Sit)]"  ^ 

dSit)  =  -at 


Js„ 


(12) 


I5(t)]" 


If  n  is  an  integer,  then  Eq.  12  becomes 


I  ij  =  -at . (13a) 


or 


if  /  \ 


5'"'V5  =  -at 


(13Z» 


For  m  7^  0,  if  #  0.  and  n  ¥=  m  (/?  and  m  are  integers).  Eqs.  13a  and  b  can 
be  integrated  as 

{[S(f)]'”'""' - 

.  lio  ' 

1 , 


E  (-1)1"  « 


/=» 

/^m—  1 


+  (-1)" 


m-1 


•Jo 


j  -  m  +  \ 
'Sit) 


In 


L^o 


=  -at. 


(14) 


If  m  is  not  an  integer,  the  solution  of  Eqs.  13a  and  b  is 


E(-i) 

/=  I 


^if- 

Oo 


{15(f)l^'"'" 
j  -  m 


-  sr"''} 

+ 1 


—at 


(15) 


The  generalized  infiltration  model  is  given  in  parametric  form  by  Eqs.  14 
and  11  or  by  Eqs.  15  and  11  for  n  being  an  integer.  For  the  cases  0  ^  n 
£  2  and  0  £  fw  £  2.  solutions  of  Eqs.  Fl  and  12  leading  to  some  special 
cases  of  the  generalized  model  are  shown  in  Table  3.  If  m  =  0  and  n  >  0. 
then  the  generalized  infiltration  model  reduces  to  the  Kostiakov  model;  and 
when  H  =  0  and  m  >  0.  the  generalized  model  becomes  the  Holtan  model, 
as  shown  in  Table  2.  If  fi  is  not  an  integer  and  m  0.  then  Eq.  12  may 
not  have  a  general  analytical  solution.  If  the  five  parameters,  a.  m.  n.  Sn. 
and  f  are  known,  then  the  infiltration  rate  can  be  calculated  numerically  by 
using  Eqs.  11  and  12  and  5  as  a  parameter  (0  £  5  £  Sn).  A  simple  least- 
squares  procedure  is  presented  in  Appendix  1  to  determine  the  five  param¬ 
eters  using  observed  data. 

Eq.  12  can  also  be  solved  in  terms  of  the  beta  functions.  Define  r  -  Sit)/ 
S„  and  Sifdr  =  dSit).  Therefore.  Eq.  12  can  be  expressed  as 


s:;- 


rfr-'dr  =  -at 


(16(7) 


or 


r  '  * 

r 

- 

s:r'"*' 

(  1  -  / 

-  Jm 

1 

1) 

-  rfdr 

( 1  bb) 
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w  >  I 


n  >  0 
n  ¥=  2 
0 
I 

■> 


Si,  -  {3ai)‘^' 

S„  -  [(/I  +  Dfl/)' 

S„e''" 

In  (5/5„)  -  (5/5„)  =  1  -uir/S,,) 

In  (5/S„)  -  2(5/5.,)  -  (l/2)(5/5„)- 
=  -((3/2)  +  (fl//5,'i)| 

[( I/5,i)  1-  at]'' 

In  (5/5„)  +  (5„/5)  =  1  -  fl, 

(5/5„)  -  (5„/5)  -  2  In  (5/5„) 

=  -(fl//5„) 

[5,',  ( 1  -  m)ctr]' 


^Ifun  equal  sign  appears  in  an  expres.sion.  that 


fc  +  a 

f.  +  ta/2)''-r'''- 
f  +  (a/9)"'/--'' 

'/  {((a  D/a'  •']/} - - 

/  Siiue'"' 

f  ^  {fa5(/)]/[5„  -  5(/)I} 

/  {a5(/)/(5„  -  5(/)]-} 

'/  *  (Va/(I/5„  ^  a/))- 
'/+  (fd  -  a/)5(/)-]/[5„  -  5(/)l} 
'/  {(a/5-'(/)l/((5„  -  5)']} 

!/  -  a(5,',""  -  ( 1  -  m)at]" 


or 


s.r 


j  (I  -  r)"'- dr  -  jT 


/'--'(I  _ 


=  at 


(16c) 


tfl ±  — fff  I  •  tinrl  #1  ““  IT  -1-  I  'T'l.  ■  I 

I.  ana  -  „  +  j  j^is  can  be  expressed  as 

s;:'"”'  ifl(w,,„*)  -  . 


(16//) 


te™7,ff';h,'  r  ‘‘1“  be  expressed  m 

terms  of  the  gamma  function  as  ^ 


fl(a;*./i*)  = 


r(//7*)r(/i*) 


H/n*  -  lit) . . 

and  BA-.)  =  tlje  incomplete  beta  function,  which  can  be  expressed  as 


Ba 


r 

a:b)  - 

Jo 


r“''(l  -  xf-'dx: 


a  >  0.  b  >  0.  r  E  (0.  1 ) 


Bjcr.b)  = 


r“(  1  -  r)" 


•-2 


-  I:y  I ) 


( 1 8c/ ) 


(18/;) 


y=()  B(a  ^  b\  j  ^  I ) 

Considerable  simplification  occurs  if ///  =  in  Eqs.  I  I  and  12.  Usine  Eq 
I.  5(D  can  be  expressed  as  an  explicit  function  of  f/o  .c 

r  •  .  "1 1  ■ 

f/'d  -/ 


5„ 


Sit)  = 


explicit  function  of  fit)  as 


(  19) 
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TABLE  4.  Special  Cases  of  the  Generalized  Model 


Parameters  in  Eq.  8 

Resulting  model 
(4) 

/n 

(1) 

It 

(2) 

f 

(3) 

0 

1 

f 

Philip  two-term  or  Green-Ampt  model 

0 

>-l 

0 

Kostiakov  model 

0 

>-l 

f 

Modified  Kostiakov  model 

I 

0 

/ 

Horton  model 

2 

0 

f 

Overton  model 

>0 

0 

f 

Holtan  model 

5=2 

Special  Cases 

Seven  of  the  eight  cited  models  are  special  cases  of  the  generalized  model. 
Some  ot  these  special  cases  can  be  inferred  from  Table  2.  For  clarity,  these 
and  other  special  cases  are  shown  in  Table  4.  The  Zhao  model  ( 1981 )  is  not 
one  ot  the  special  cases  for  the  reasons  cited  earlier. 

For  m  =  0.  n  >  0.  Eq.  12  becomes 

S(i)  =  S„  -  |(/i  -  Dar)'  """  . 

Using  Eqs.  8  and  9.  Eq.  20  leads  to 


Eq.  21  is  the  modified  Kostiakov  model.  If/  =  0.  Eq.  21  reduces  to  the 
Kostiakov  model. 

Similarly,  for  w  =  0.  /t  =  I,  the  generalized  model  reduces  to  the  Philip 
two-term  model  (Philip  1957)  or  the  Green-Ampt  model  (Green  and  Ampt 
1911).  and  also  to  the  Kirkham-Feng  model  (1949)  for  horizontal  How.  It 
will  be  shown  in  the  ensuing  discussion  that  the  Green-Ampt  model  is  equiv¬ 
alent  to  the  Philip  two-term  model.  If ;?  =  0.  m  =  1.  the  aeneralized  model 
becomes  the  Horton  model  (1938).  and  for  /i  =  0.  w  =  2.  the  resuitina 
model  is  ttye  Ovenon  model  (1964). 

For  m  >  L  /?  =  0.  and  /n  ^  2.  Eq.  12  has  the  solution 

S(r)  =  -r  (m  -  Dar)"'’”"" . ^22) 

Thcrctore.  Eq.  8  becomes 

An  cjS^(f) . 

Substituting  Eq.  22  m  Eq.  23  yields 

/(n  =y;.  -  -  (1  -  . (24) 

Eq.  24.  the  general  solution  for  the  Holtan  model  (1961).  docs  not  appear 
to  have  been  reported  before. 

.S(’mc  ot  the  intiitration  models  listed  in  column  5  of  Table  2  may  not 
lo('K  e.xactly  like  their  usually  reponed  form  shown  in  column  6  of  the  same 
tabic.  The.se  torms  ot  infiltration  models  arc  actually  equivalent.  For  e.\- 
amplc.  the  Horton  model  i  1938)  in  the  current  tbrm  is 


(20) 

(21) 
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(25) 


5o 

fit)  =f^  +  -e-"‘‘ . 

a 

If  we  let  /(O)  =  /o.  Eq.  25  yields 
So 

-=fo-f . (26) 

a 

Substitution  of  Eq.  26  into  Eq.  25  produces  exactly  the  same  form  of  the 
Horton  model. 

fU)  =  f  +  (fo  . (27) 

Some  of  the  models  in  Table  2  do  not  have  general  analytical  solutions  (Hol- 
tan  model)  or  do  not  have  explicit  analytical  solutions  for /(/)  (Green- Ampt 
model)  as  reported  in  literature.  However,  simpler  and  explicit  solutions  of 
fU)  can  be  obtained  using  the  systems  approach  for  all  of  the  eight  models. 


Equivalence  Between  Green-Ampt  and  Philip  Two-Term  Models 

The  Green-Ampt  and  the  Philip  two-term  models  have  the  same  form  of 
solution  in  column  5  of  Table  2.  Since  h,  9,.  and  5,  are  all  constants  in  the 
Green-Ampt  model,  if  K  =  f  and  a  =  KB,(h  +  5,  ).  the  Green-Ampt  and 
the  Philip  two-term  models  are  equivalent.  From  Table  2.  these  two  models 
have  the  same  expression  for  SO)' 

SU)  =  S„  -  VZar . (28) 

Substituting  Eq.  28  into  Eq.  2.  and  solving  with  /,  =  /  =  AT.  we  have 


/(f)  =  /  +  a 


(29) 


which  is  the  Philip  two-term  infiltration  model. 

On  rhe  other  hand.  Eq.  28  can  be  expressed  in  the  form  of  the  Green- 
Ampt  model.  We  rewrite  Eq.  28  as 

l5o  -  5(f)]'  =  2ar . (30) 

Eq.  30  can  also  be  expressed  as 

. (31) 

Jo  Jo 


or 

|5(,  -  5(f)]c/15o  -  Sir)}  =  adt  . (32) 

Rearranging  Eq.  32  yields 
d\Sn  -  5(r)|  a 

-  =  .  . (33) 

dr  5f,  -  5(/) 

Substituting  Eq.  33  into  Eq.  2,  replacing  /'  by  K  and  solving,  we  obtain 


fir) 


(34) 
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Notice  that  So  -  Sit)  =  W(t)  =  L6,,  where  L  is  the  depth  of  the  transmission 
zone,  and  0,  is  volumetric  moisture  content.  Eq  34  can  then  be  written  with 
a  =  K%,ih  +  S<  )  as 


fit)  =  K 


1  + 


9,(/t  +  S,) 

W 


=  K 


L  +  h  + 


(35) 


which  is  the  Green- Ampt  model. 

From  this  proof,  we  see  that  the  hydraulic  conductivity  K  defined  by  Green 
and  Ampt  (1911)  has  the  same  value  as  the  steady-state  infiltration  rate  in 
the  Philip  two-term  model  (Philip  1957).  The  coefficient  a  in  the  Philip  two- 
term  model  is  proportional  to  the  sum  of  the  hydraulic  head  h  and  the  cap¬ 
illary  suction  head  S,  .  Again,  the  solutions  by  systems  approach  are  of  sim¬ 
pler  form  than  those  derived  using  some  of  the  other  methods. 


Special  Case  of  General  Model 


It  may  be  useful  to  e.xamine  a  special  case  of  the  general  model.  For  the 
case  n  =  1.  w  =  1.  the  generalized  model  (Eq.  11)  becomes 


fit)  =/  + 


aSit) 

Sn  ~  Sit) 


(36) 


and  Sit)  is  determined  by  integration  of  Eq.  12: 


t 


So  S 
^  .So 


(37) 


Eqs.  36  and  37  can  be  considered  as  a  parametric  model  of  fit)  with  S  as 
the  characteristic  parameter.  This  model  is  actually  a  combination  of  the 
Horton  model  and  the  Green- Ampt  model  (or  the  Philip  two-term  model). 
When  Sit)  is  removed  from  the  numerator  of  the  second  term  on  the  right- 
hand  side  of  Eq.  36.  this  equation  becomes  the  Philip  two-term  model;  and 
when  the  term  So  —  S(r)  is  removed  from  the  denominator.  Eq.  36  becomes 
the  Horton  model. 

The  parametric  model  of  Eqs.  36  and  37  can  also  be  expressed  in  implicit 
form  with  i*  as  a  function  of  fit).  From  Eq.  36 


S(/)  = 


a  ^  fit)  -f. 


(38) 


Substitution  of  Eq.  38  in  Eq.  37  leads  to 

f  fit)-f  .  r  f(t)-f. 

r  =  —  s - In  - 

a  U  fit)  -  f  -  fit)  -  f  _ 

This  model  consists  of  only  three  parameters;  f  .  S„.  and  a.  It  is  not  only  a 
simple  infiltration  model,  but  also  includes  consideration  of  the  soil-water 
content  and  movement. 

Following  the  concept  embodied  in  this  infiltration  model,  it  may  be  pos¬ 
sible  to  combine  existing  simple  models  to  derive  from  the  general  model 
new  infiltration  models  with  three  or  less  parameters.  The  parameters/  and 
5„  of  the  new  models  can  usually  be  obtained  from  parameters  evaluated  for 


847 


TABLE  5.  Relationship  among  Parameters  in  Derived  and  Reported  Forms  of 
Infiltration  Models 


Model 

(1) 

Parameters  in 
reported  form 
(2) 

Parameters  in 
current  form 
(3) 

Relationship 

(4) 

Honon 

/o 

S{) 

So  -  at  /,  -  /  ) 

Honon 

L 

f. 

/  =/ 

Honon 

a 

a 

a  -  a 

Ovenon 

f- 

/ 

f.  =f 

Ovenon 

a 

a 

a  =  a 

Ovenon 

5„ 

s„ 

S{)  “ 

Philip 

A 

f 

II 

Philip 

s. 

a 

a  =  (5:)/2 

Kostiakov 

K 

n 

n  =  A/{\  -  A) 

Kostiakov 

A 

a 

a  =  AK'  ' 

Modified  Kostiakov 

f 

f 

/  =/ 

Modified  Kostiakov 

K 

n 

/;  =  A/{  1  -  A) 

Modified  Kostiakov 

A 

a 

w  =  A/r'  ^ 

Zhao 

a 

a 

a  =  a 

Zhao 

b 

b 

b  b 

Zhao 

Wo 

— 

— 

other  models  as  shown  in  Table  5  or  estimated  directly  from  observ'ed  data. 
The  proportionality  factor  a  may  be  directly  determined  from  observ  ed  data 
using  the  method  of  least  squares  or  moments. 

First-Order  Analysis  of  Uncertainty 

Eq.  8  is  the  fundamental  equation  of  the  generalized  infiltration  model 
proposed  here.  There  may,  however,  be  uncertainties  involved  with  this 
equation  that  must  be  recognized.  Besides  the  uncertainty  of  the  equation 
itself,  there  can  be  two  other  sources  of  error  in  / due  to  spatial  and  temporal 
variability.  First,  the  measurement  of  S  may  be  inaccurate.  Second,  the  pa¬ 
rameters,/  'n,  and  n  may  be  subject  to  errors.  The  first-order  analysis 

of  uncertainty  can  be  employed  to  quantify  the  expected  variability  arising 
from  uncertainty  in  5  and/or  a,/f..  So.  'n,  and  n.  To  illustrate  the  procedure, 
let  us  represent  Eq.  8  as 

f  =  g{S,a,f,..So.m.n) . (40) 

where  g  =  some  function.  It  is  implied  in  Eq.  40  that  /and  5  depend  on  /. 
If  the  true  values  of  S,  a,/.,  and  //  differ  from  their  corresponding 


nominal  (average)  values.  5.  d,  /  ,  So.  m.  and  n,  where 

/=  g(S.d./ ,So./^.//) . (41) 

then  the  effect  of  this  discrepancy  on  /  can  be  evaluated  by  expanding  the 
function  g{  •  )  around  the  point  defined  by  S  =  S.  u  =  d.  i  =  /  .  S.  =  S*. 
m  =  m.  n  =  n. 

For  purposes  of  simplicity,  let  .r,  =  S.  .V:  =  a,  f  .  -v,  =  S,.  -x\  =  m. 
and  =  n.  Then  the  Taylor  series  expansion  of  Eq.  40  yields 
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/  =  g{Ld.f,.J„,m.n) 


(-r,  -  X.) 


+ 


(.V,  —  .v,)(.v,  -  Xj)  +  higher-order  terms 


(42) 


where  [(dg/3,r,)-]  =  the  first-order  partial  derivative  of  the  function  g  with 
respect  to.r,  evaluated  at  the  point  of  expansion  and  (d\g/ 

5.t,d.r, )-  =  the  second-order  partial  derivative  of  g  with  respect  to  v,  and  .v, 
evaluated  at  the  point  of  expansion. 

If  the  second-  and  higher-order  terms  are  neglected,  the  resulting  first- 
order  expression  for  the  error  in  /  is 

/-/=£  =  2  (^)  (-t,  -.f,) . (43) 


The  variance  of  this  error  is 

S;  =  E[{f-ff\ . (44u) 


i44b) 


where  E  =  the  expectation  operator.  For  simplicity,  it  is  assumed  that  the 
variables  x,  are  independent. 


or 


where  Si  =  the  variance  of  .v, .  From  Eq.  8.  the  partial  derivatives  ic»e/ii.v, 
=  df/dx,)  arc 


d/  amS"  '  anS"' 

BS  (So  -  5)"  (So  -  S)"* ' 


( 46n  I 


^  -  ■S'" 

oa  (S)  ~  .S)" 

—  - 

Bf  anS”' 

BSo  (S,  -5)'"' 


(46/7) 


(46c) 


(46c/) 
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(46e) 


1. 

dm 

dn 


=  (/-/c)  log  S . 

=  -if-fc)  log(5o  -  5) 


(46/) 


where  log  denotes  the  natural  logarithm. 

By  substituting  the  values  of  these  partial  derivatives  in  Eqs.  45a  and  b, 
the  error  in  /  can  be  computed.  If  the  error  in  /  is  computed  due  to  error  in 
only  one  variable,  say  5,  then 


S}  =  var  (/)  =  (  ^  )  S|  = 
do 


amS”~' 


1  + 


nS 


(5o  -  Sr  L  m(So  -  S)J 


(47) 


The  value  of  Sj  =  the  standard  error  of  estimate  of  /.  If  error  in  /is  calculated 
due  to  error  in  a  only,  then 


var  (/)  =  S;  =  (  ^  )  5;  = 
da 


(So  -  srj 


s: 


If  only  /.  is  to  be  considered,  then 


. 

If  the  variable  under  consideration  is  So,  then 


Sf  = 


dSoJi, 


anS” 


L(So  -  S)’'"'J 

If  the  variable  under  consideration  is  m.  then 


(48) 


(49) 


(50) 


5;  =  (  —  S-  =  [(/-/.)  log  Sy-Sl . .  .  (51) 

\dm/ 


If  the  variable  under  consideration  is  n,  then 


^g 


S;  =  —  5;  =  [-(/  - /)  log  (5o  -  S)]=S; 

\dn/ 


(52) 


The  total  variance  of  infiltration  rate  can  also  be  expressed  as  a  function 


var  (f)  =  Sf  = 


amS" 


(So  -  sr 


1  + 


nS 


m(S„  -  S)J 


So,  m. 

and  n  as 

11 

U  - 

5" 

Jj 

I  ’ 

.(So  -  sr. 

S: 


^  ^  +  [(/-/  )  log  S]^5; 


^  [-(/-/.)  log(5o  -  S)]‘S; . (53) 

The  variance  of  infiltration  rate  can  also  be  expressed  as  a  function  of  the 
coefficients  of  variation  CV  of  5,  a,/,  5o,  m.  and  n  as 


var  ( /)  =  f'CV'i  f)  = 


ay" 


(5o  -  sr] 


m 


1  - 


nS 


mi  Si)  -  5)J 


-  CV: 
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F!w  3.  Standard  Error  of  Estimate  oi  fu)  as  Function  of  Reiative  Deviations  in 
Parameter:^  =  0.1.  ,v  =  1,  =  2.  5  =  0.3,  So  =  0.6,/  =  0.001;  1  =  Deviation  in 

Variable  5  Only;  2  =  Deviation  in  Parameter  A  Only;  3  =  Variation  in  Parameter 
FC  Only;  4  =  Deviation  in  Parameter  S„  Only;  5  =  Deviation  in  Parameter  A/  Only; 
and  6  =  Deviation  in  Parameter  .v  Only 


^  (n5„)-CVi,  -  (w  log  srcvi  *  [/,  log  (So  -  5)]-CV,y  -  ffCVj . (54) 

It  is  seen  from  Eq.  53  that  the  variance  of  fit)  not  only  depends  on  the 
variance  of  each  parameter  but  also  on  the  values  of  the  parameters  them- 
^lves._  An* illustrative  example  is  shown  in  Fig.  3  for  d  =  0.1.  «  =  \.  m 
~  2,  S  —  0.3  m.  5,1  =  0.6  m.  and/,.  =  0.001  m/h.  These  values  represent 
normal  cases.  Fig.  3  shows  that  the  standard  error  of  estimate  of  fit)  is  very 
sensitive  to  the  reiative  deviations  in  parameters  5.  m.  5o,  a.  and  n  and  is 
relatively  insensitive  to  the  relative  deviation  in  parameter  /  .  This  conclu¬ 
sion.  however,  is  based  on  the  mean  values  used  and  mav  chanee  with  the 
changes  in  the  mean  values. 


Verification  of  Model 

The  generalized  model  has  five  parameters;  /  (steady-state  infiltration  rate). 
5„  (initial  storage  space),  a  (a  proportionality  factor),  and  m  and  n  (two 
exponents).  Table  5  presents  the  relationships  between  the  parameters  of  the 
infiltration  models  derived  by  the  systems  approach  listed  in  column  5  of 
Table  2  and  those  reponed  in  literature  for  all  six  explicit  and  analytical 
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FIG.  4.  Calibration  of  Infiltration  Model:  Soil  Type:  Alapaha  Loamy  Sand  ID  = 
01 011 D;  and  Soil  Type:  Cowarts  Loamy  Sand  ID  =  17032D  (>*'  -  Observed  Data 
Set  1 ;  ^  =  Observed  Data  Set  2;  and - -  Calculated  Value) 

infiltration  models  listed  in  column  6.  Once  the  parameters  of  the  usual  form 
of  a  listed  model  are  recognized,  the  corresponding  parameters  for  the  model 
derived  by  the  systems  approach  can  be  easily  found. 

If  infiltration  observations  are  available  for  a  soil,  then  the  five  model 
parameters  can  be  estimated  for  that  soil  using  an  iterative  procedure.  How¬ 
ever,  this  procedure  is  usually  complicated  and  does  not  take  into  account 
the  physical  import  of  the  parameters.  A  simplified  procedure  was  therefore 
employed  in  this  study.  The  procedure  is  based  on  the  premise  that/  and 
So  can  be  initially  determined  from  the  knowledge  of  soil  propenies,  ante¬ 
cedent  soil  moisture  conditions,  and  infiltration  measurements.  This  premise 
is  not  unduly  restrictive,  for  both  /  ^^nd  So  are  physical  quantities.  On  one 
hand,  /.  is  the  ultimate  infiltration  rate  and  can  be  determined  experimen¬ 
tally.  In  practice./,  should  satisfy 

0  </  <  min  l/o(r)]  . (55) 

where /,(r,)  =  the  iih  obser\’ed  infiltration  rate:  and  5,)  =  the  effective  po¬ 
tential  storage.  Experimentally,  it  is  equal  to  the  soil  column  times  the  ef¬ 
fective  porosity.  From  Eq.  II.  if  the  soil  column  is  saturated  5  =  0.  then 
infiltration  rate  will  equal  /  .  At  /  =  0,  5  =  5n.  then  /(O)  —  This  is 
consistent  with  the  Horton  model.  In  practice.  5n  can  be  selected  by  con- 
sidenng  a  practical  soil  depth  and  an  average  soil  porosity,  since  only  the 
top  soil  zone  may  significantly  affect  the  dynamic  infiltration  process  (Fok 
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RG.  5.  Calibration  of  Infiltration  Model:  Soil  Type:  Stilson  Loamy  Sand  ID  =  101 01 D 
=  Observed;  -  =  Calculated;  =  0.00093;  B  =  1.28:  V  =  0.353;  5..  =  27.43 
(cm);  and  FC  =  3.048  (cm/h)] 


and  Chiang  1984).  Normally,  this  top  soil  zone  is  less  than  2  m.  If  average 
soil  porosity  is  0.3.  then  5,,  should  satisfy 

0  <  5o  <  0.6  . . . (56) 

Ovenon  (1964)  assumed  =  0,2743  m  (soil  depth  =  36  in.),  while  Holtan 
( 1961)  chose  5,,  =  0.0457  m  (soil  depth  =  6  in.).  On  the  other  hand,  it  may 
be  difficult  to  estimate  accurately  the  actual  values  of  these  two  parameters 
in  a  watershed  due  to  spatial  and  temporal  variabilities.  However,  once  / 
and  So  are  reasonably  chosen  following  Eqs.  55  and  56.  then  the  other  three 
parameters  a.  m,  and  n  can  be  estimated  using  the  Icast-.squares  method. 
Hence,  a  good  tit  to  tield  data  can  be  achieved.  An  explicit  least-squares 
solution  for  estimation  ot  these  parameters  is  given  m  Appendix  1. 

Ten  obser\'cd  infiltration  data  sets  from  Rawls  et  al.  (1976)  were  used  to 
test  the  goodness  ot  the  generalized  infiltration  model  given  m  this  paper. 
The  generalized  model  usually  gave  a  belter  fit  than  other  special  models 
tor  each  data  set.  The  computational  procedure  Is  given  in  Appendix  1.  For 
three  sample  data  sets,  observed  and  computed  infiltration  curves  arc  given 
in  Figs.  4  and  5.  respectively,  for  different  types  of  soil. 
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Because  the  five  parameters  can  be  determined  directly  from  Appendix  I, 
no  iteration  or  trial-and-error  method  is  involved.  However,  optimal  sets  of 
parameters  may  not  always  be  unique,  and  more  than  one  set  of  optimal 
parameters  may  lead  to  an  equally  good  fit  to  observed  data.  From  physical 
considerations,  it  may,  therefore,  be  preferable  to  choose  0.25  m  :£  5o  ^ 
0.45  m,  and as  slightly  smaller  than  the  smallest  observed  infiltration  rate 
value.  Then  parameters  a,  m,  and  n  are  determined  by  Appendix  I.  If  the 
computed  value  of  a  is  too  small  (o  <  10“’  in  cm-h  units),  one  can  change 
fc  slightly  in  the  range  given  by  Eq.  56.  Normally,  one  or  two  trials  will 
produce  a  reasonable  set  of  optimal  parameters.  In  the  three  data  sets  shown 
in  Figs.  4  and  5,  5©  was  fixed  at  27.43  cm  where  the  soil  depth  considered 
is  91.44  cm  and  soil  porosity  is  30%.  In  the  range  given  in  Appendix  I,/, 
was  selected  so  that  a  was  in  the  range  from  10“’  to  lO"**.  All  computations 
showed  that  the  five  parameters  were  easily  determined,  and  hence  produced 
a  good  infiltration  model. 

Concluding  Remarks 

This  study  utilizes  systems  approach  to  derive  a  generalized  infiltration 
model.  Seven  popular  infiltration  models  are  shown  to  be  special  cases  of 
this  generalized  model.  This  model  has  five  parameters;  steady-state  infil¬ 
tration  rate  initial  storage  space  available  So;  a  proportionality  factor  a; 
and  two  exponents  m  and  n.  It  may  be  useful  to  combine  two  or  three  pop¬ 
ular  infiltration  models  into  one,  according  to  the  generalized  model,  and 
yield  a  three-parameter  model  as  illustrated  by  Eqs.  36  and  37.  For  10  ob¬ 
served  data  sets,  the  procedure  for  parameter  estimation  was  quick  and  ef¬ 
ficient.  The  model  reproduced  infiltration  rates  for  these  data  sets  reasonably 
well. 

Acknowledgments 

The  work  reported  in  this  paper  was  supported  in  part  by  funds  provided 
by  the  Army  Research  Office  through  Southern  University  under  the  project 
on  “A  Continuum  Model  for  Streamflow  Synthesis." 

I 

Appendix  I.  Least-Squares  Estimation  of  Model  Parameters 

It  is  assumed  that  the  model  parameters  and  So  are  obtained  from  soil 
properties  and  antecedent  conditions.  The  parameters  estimated  are  a.  m. 
and  n.  To  that  end.  the  least-squares  objective  function  is  constructed  in  log 
space  as 

V 

55  =  2  {log  l/cbs(/)  -  log  UoJO  -/.]}’ . (57) 

/=  I 

where  /ohs(0  =  observ^ed  infiltration  rate  at  the  /th  time;  =  computed 

infiltration  rate  at  the  /th  time:  55  =  the  error:  and  /V  =  the  number  of 
observations  or  times.  It  is  recognized  that  the  sum  of  squares  of  log  values 
is  minimized,  not  the  sum  of  the  infiltration  rates  themselves.  As  a  result, 
the  estimate  of  a,  n,  and  m  may  be  somewhat  biased.  Minimization  based 
on  the  log  values  was  preferred  for  two  reasons.  First,  it  leads  to  an  ana- 
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lytical  solution.  Second,  experimental  verification  of  the  model  showed  that 
the  model  fit  the  field  data  reasonably  well.  Minimization  based  on  arith¬ 
metic  values  of  infiltration  rates  leads  to  equations  that  are  not  explicit  in 
terms  of  a,  n,  and  nt,  and  hence  an  iterative  solution,  although  not  very 
difficult  to  accomplish  with  a  computer,  is  required. 

Eq.  57  can  be  written  as 

=  i  (log  eg  . ,58, 

where  S(/)  =  the  /th  value  of  S(t).  For  simplicity  of  notations,  the  suffix  / 
is  dropped  from  now  on,  keeping  in  mind  that/ob,  and  5  both  are  functions 
of  time  and  vary  from  one  observation  to  the  other.  Eq.  58  is  written  as 

55  =  2  [log  (/obs  -/.)  -  log  a  -  m  log  5  +  n  log  (5o  -  5)f . (59) 

The  parameters  a.  m,  and  n  are  obtained  so  as  to  produce  the  minimum 
value  of  55.  To  that  end,  55  is  partially  differentiated  with  respect  to  each 
parameter  separately,  and  each  partial  derivative  is  equated  to  zero.  This 


yields 

355  1 

-  =  0  =  —  [2  log  (f^  -/,)  -N\oga-mL  log  5 

da  a 

+  log  (5o  -  5)] . (60) 

or 

Af  log  a  +  m2  log  S  -  ril  log  (5o  -  5)  =  2  log  -  X) . (61a) 

355 

T  “0 - ^  log  (/obs  ~  fc)  log  5  +  log  al  log  5  +  m2  (log  5)‘ 

dm 

-  /i2  log  (5o  -  5)  log  5 . {6\b) 

or 

log  a1  log  5  +  m2  (log  5)*  -  n1  log  (5o  -  5)  log  5 

=  2  log  (/obs,-/)  log  5 . (62a) 

355 

-  0  -  2  log  (/obs  -  fc)  log  (5o  -  5)  -  log  a2  log  (5o  -  5) 
dn 

-  m2  log  (5o  -  5)  log  5  n2  [log  (5o  -  5)]’ . (626) 

or 

log  a2  log  (5o  -  5)  +  m2  log  (5o  -  5)  log  5 

-  nl  [log  (5„  -  5)]^  =  2  log  (/obs  -/.)  log  (5o  -  5)  . (63) 


Eqs.  61-63  are  solved  simultaneously  by  Kramer's  rule  to  obtain  the  pa¬ 
rameters  a.  m,  and  n.  In  matrix  form,  these  equations  can  be  written  as 


N 

2  log  5 

-2  log  (5o  -  5) 

log  a~ 

2  log  5 

2  (log  5)' 

-2  log  (5o  -  5)  log  5 

m 

2  log  (5o  -  5) 

2  [log  (5o  -5)1  log  5 

-2  [log  (5o  -  5)/ 

n 
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(64) 


s  log  {fobs  - fc) 

=  2  log  (/ob3  - /J  log  5 

_1  log  (/obs  -  fc)  log  (5o  -  S) 

For  notational  convenience  and  writing,  let  /I  =  2  log  5,  5  =  2  (log  5)% 
C  =  2  log  (So  -  5),  D  =  2  [log  (5o  -  5)]  log  5,  G  =  2  [log  (5o  -  5)]'. 

-  2  log  (/obj  -  /  =  2  log  (/ob,  -  /.)  log  5;  and  y  =  2  log  ( /.b,  - 

/ )  log  (5o  —  5).  Eq.  64  can  be  rewritten  using  these  notations  as 


~N  A  -C~ 

log  a 

A  B  -D 

m 

= 

1 

C  D  -G 

n 

J 

Then 


1 

log  a  =  — 

D. 


H  A  -C' 
I  B  -D 
J  D  -G 


where  D,  —  the  determinant  of  the  coefficient  matrix  expressed  as 

D.  = 


N  A  -C 
A  B  -D 
C  D  -G 


D,  =  N(D-  -  BG)  -  A(CD  -  AG)  -  C(AD  -  CB) 
Hence 

H(D-  -  BG)  -  AUD  -  IG)  -  CUD  -  JB) 


loe  a  = 


m 


n  = 


N  H  -C 
A  I  -D 
C  J  -G 

N  A  H 
A  B  I 
C  F  J 


D, 


NUD  -  IG)  -  H(CD  -  AG)  -  C(Ay  -  Cl) 

D, 


N(BJ  -  DI)  -  A(AJ  -  Cl)  ^  HiAD  -  CB) 


D, 


(66) 


(67a) 

(676) 

.  (68) 

-  (69) 

(70) 


Eqs.  68-70  explicitly  yield  optimal  estimates  of  a.  m.  and  n. 

For  a  given  set  of  infiltration  data,  i.e.,  [f,,/obs(/)].  /  =  1.2 . N.  one 

can  determine  the  five  parameters  as  follows. 


Step  1 

Assume  initial  values  r  =  0.  5(0)  =  5n. /(O)  =/ob,(l). 

Step  2 

Based  on  the  physical  meaning  of  parameters  5„  and/  and  observed  data, 
choose/  in  the  range 

0  </  <  min  l/,bs(/)l.  /  =  1.2 . N . (71) 

where  /,bv(/)  =  the  obser\ed  infiltration  rate.  For  most  types  of  soils  the 
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infiltration  rate  may  be  significantly  affected  only  by  the  top  soil  zone  (say. 
soil  depth  less  than  2  m).  If  soil  porosity  is  30%.  then  So  may  be  chosen  in 
the  range 

0  <  5o  <  0.6  m . (72) 


Step  3 

From  Eq.  2  or  Eq.  1.  one  can  calculate  S(r)  as 

s,..i  =  s,  +  {/  -  o.5[y„bs(/  +  1)  +/o( /)]}(?,... I  -  f,) 
s,  =  5(0)  =  5...  /  =  1 .  2 . N . 


Step  4 

Compute  parameters  a,  m,  and  n  by  Eqs.  68.  69.  and  70.  respectively. 


Step  5 

Once  all  five  parameters  are  determined,  the  values  of  infiltration  rate  can 
be  generated  by  the  model  for  a  given  set  of  values  of  5, .  0  s  5,  <  5„.  Thus 


f.  =/;.  - 


aS'" 

(So  -  S, )" 


(74) 


and 


(5,.  I  -  5,) 

/  -  o.5(,/;.,  ^/y 


i  =  1.2 . ;V 


(75) 
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Appendix  III.  Notation 

The  following  symbols  are  used  in  this  paper: 

a  =  proportional  factor  of  infiltration  model; 

F  =  cumulative  infiltrated  water  depth  (cm,  L); 

/  =  infiltration  rate  (cm/h,  L/T); 
f  =  steady-state  infiltration  rate  (cm/h,  L/T); 
f  =  seepage  rate  (cm/h,  L/T); 

h  =  ponding  water  depth  (cm,  L); 

K  =  hydraulic  conductivity  in  transmission  zone  (cm/h,  L/T); 

L  =  soil  depth  (cm.  L); 

m  =  exponent  (dimensionless); 

n  =  exponent  (dimensionless); 

S  =  potential  storage  (cm  or  cm\  L  or  L’); 

Sr  =  capillary  suction  head  (cm,  L); 

5o  =  initial  soil  space  available  (cm  or  cm\  L  or  L^); 

t  =  time  (min.  T); 

W  =  cumulative  infiltrated  water  depth  (cm.  L);  and 
9,  =  volumetric  soil  water  content  (cmVcm\  LVL')- 
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